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ABSTRACT 


State Estimation of International Space Station Centrifuge Rotor with Incomplete 

Knowledge of Disturbance Inputs 

by 

Michael James Sullivan 

This thesis develops a state estimation algorithm for the Centrifuge Rotor (CR) system 
where only relative measurements are available with limited knowledge of both rotor 
imbalance disturbances and International Space Station (ISS) thruster disturbances. A 
Kalman filter is applied to a plant model augmented with sinusoidal disturbance states used 
to model both the effect of the rotor imbalance and the ISS thrusters on the CR relative 
motion measurement. The sinusoidal disturbance states compensate for the lack of the 
availability of plant inputs for use in the Kalman filter. Testing confirms that complete 
disturbance modeling is necessary to ensure reliable estimation. Further testing goes on to 
show that increased estimator operational bandwidth can be achieved through the 
expansion of the disturbance model within the filter dynamics. In addition, Monte Carlo 
analysis shows the varying levels of robustness against defined plant/filter uncertainty 


variations. 
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1 Introduction 

The National Aeronautics and Space Administration (NASA), in association with the 
Japanese Aerospace Exploration Agency (JAXA), are building a Centrifuge 
Accommodation Module (CAM) for attachment onto the International Space Station (ISS). 
The CAM houses the Centrifuge Rotor (CR) and will be attached at node 2 on the 
International Space Station (ISS) as shown in Figure 1-1. 



Figure 1-1. Exploded Diagram of International Space Station Components [1] 
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The CAM/CR is an orbiting laboratory which will study the effects of zero gravity and 
micro gravity environments on rodents. More details concerning the CAM/CR can be 


found in Section 1.1. 
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Before the CAM/CR can be attached to the ISS, all verification must be completed on the 
ground to ensure robust stability and safe operation. This is important not only to for 
increasing the probability of mission success, but also to make certain the safety of the ISS 
crew members. One issue concerning the safe operation of the CR aboard the ISS is the 
occurrence and the effect of rotor imbalances due to changing inertia during CR operation. 
The effect of rotor imbalances can be found in everyday life such as an unbalanced 
washing machine drum impacting the side of the washing machine or steering issues 
caused by unbalanced automotive tires. Although this problem may seen benign in the 
washing machine example, if the massive rotor in the CR impacts the CAM, critical 
damage to the ISS could result. 

One method of solving this problem would be to use counterbalancing masses to cancel out 
any imbalances in the rotor. In the case of the CR, a system called the Auto Balancing 
System (ABS) employs this method. A problem occurs during implementation of the ABS 
due to the unavailability of measurements integral to ABS control, namely the rotor’s 
absolute (i.e., relative to inertial space) states. These absolute rotor states cannot be 
obtained, because the displacement sensors are located in such a manner that only relative 
(i.e., between two moving masses) measurement are possible. Therefore an estimator is 
needed to estimate absolute rotor states from the relative measurements. This thesis 
proposes a method for estimating absolute rotor states from available relative/corrupt 
measurements involving the use of a Kalman filter. However using a standard Kalman 
filter formulation requires availability of both the rotor imbalance disturbances as well as 


the ISS thruster disturbances which are not available. This thesis will also discuss the 
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methods used to overcome this problem. Note that the words filter, observer, and estimator 
will be used interchangeably throughout the thesis. 

Although the goal of this thesis is specific, the basic premise of the problem being solved is 
applicable to any field where there is a need to compute absolute measurements from 
relative and/or corrupt measurements with limited input knowledge. For example, state 
estimation would be helpful in many applications such as determining the core temperature 
of a nuclear reactor, where it is too hazardous for sensor location. This is accomplished 
with the use of thermodynamics and sensors placed in less intense locations [3]. Also, 
optimal filters are useful in the field of aeronautics when applied to estimation of turbine 
blade states through dynamics and inferior measurements [4]. 


1.1 Introduction to Centrifuge Accommodation Module (CAM) 

The CAM, shown in Figure 1-2, which is composed of a life sciences glove box and freezer 
racks also houses the CR. The CR contains up to 4 habitats designed to house rodents. 
The CR will be used to study the long term effects of zero gravity and micro gravity 
environments on rodents. An artificial gravitational force of anywhere from 0 to 2 g can be 
generated by spinning the rotor anywhere from 0 and 1.4 Hz. The normal operational spin 
rate is 0.7 Hz, which is the spin rate necessary for 1 g. 
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Figure 1-2. CAM Internal Components [2] 


A rotor imbalance will occur whenever the spinning-member center of mass is not on the 
spin axis (e.g., due to location of rodents). Also a disturbance caused by the ISS jet-firing 
Attitude Control System will act on the rotor through the CAM shroud. Two separate 
systems, the Vibration Isolation Mechanism (VIM) and the Auto Balancing System (ABS) 
will be used to help minimize the rotor motion caused by these two disturbance sources. 
They are shown in Figure 1-3. Excessive rotor motion will result in snubber strikes against 
the shroud, causing the system to perform a safety shutdown. 
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Figure 1-3. CR Components 


The ABS controls the counter balancing masses which move in order to cancel out any 
rotor imbalance caused by rodent motion [5]. The sensors, which measure the motion of 
the rotor relative to the ISS, are located within the VIM. This relative measurement is the 
only available measurement with information of rotor motion. 
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By sensing only relative motion, if any ISS motion occurs, the relative measurement will 
not be the resulting motion due to pure rotor imbalance. In the case of a balanced rotor 
with only ISS motion, if the ABS were to act on the relative motion alone, it would drive 
the balancing mass away from the spin axis, effectively introducing an imbalance into a 
previously balanced rotor. This is shown in Figure 1-4. 

The CR controls and sensors do not interface with the ISS controls and sensors. This lack 
of system interaction limits the amount of knowledge available for either system’s 
controllers. The result is no direct knowledge of the ISS disturbance inputs which affect 
the relative measurement sensor located within the VIM. Also, since rodent motion is 
unpredictable and unmeasured, neither a rotor disturbance measurement, d r , nor an ISS 
disturbance measurement, d s , is available for use by the ABS controller or for use by the 
Kalman filter during state estimation. However, some rotor and ISS disturbance 
parameters (spin frequency, ISS disturbance characteristics, approximant rodent mass, etc.) 
are nominally known. Furthermore, the measurement, x rel is corrupted by the addition of 

sensor noise, v k . 



Figure 1-5. Overall CR System and Control 
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Since the main purpose of the VIM is to allow the rotor to follow the rigid body motion of 
the ISS while isolating rotor vibrations, a relative measurement is sufficient for the 
Vibration Isolation Controller (VIC) (see Figure 1-5); this can achieved by ensuring there is 
no change in relative displacement. However, the main purpose of the ABS is to 
counterbalance any imbalance caused solely by rodent motion, therefore, a relative 
measurement is not sufficient. Instead, absolute rotor state information is necessary for 
proper ABS control. This leads to the central question addressed in this thesis; that is, 
“How do we calculate absolute rotor states from relative measurements, with only partial 
knowledge of the disturbance inputs into the system?” This thesis provides a method of 
state estimation through the use of a Kalman Filter applied to a plant model which has been 
augmented by disturbance states. A more in depth discussion on the system and the details 
of the CR example problem can be found in Chapter 2. 



Figure 1-6. Open Loop System Used for Thesis 
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Since state estimation is the only operation of concern, the controllers in Figure 1-5 will be 
eliminated to create the open loop system found in Figure 1-6. This is the system used 
during the filter design process. 

1.2 Alternative Estimation Options 

One approach to resolve the lack of input knowledge is input reconstruction. Input 
reconstruction involves the use of the knowledge of the plant and the output time history to 
estimate the input, u, which in this case would include both rotor and ISS disturbances. 
This is also known as Inverse System Identification technique. This method may be 
helpful during state estimation, because if the inputs into the system can be reconstructed, 
then they can be used in the estimation process (see Figure 1-7). 



Figure 1-7. Possible Use of Inverse System ID for Estimation Process 
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A time domain method for estimating the applied forces on a structure was proposed by 
Stelzner, Kammer, and Milenkovic [6]. The method uses a non-causal moving average 
representation of the inverse structural system and has been successful in estimating the 
individual input forces for structures where the sensors are not collocated with the force 
input locations. The problem with the implementing of this method is that it only allows 
for near - real time estimation of the input forces, while the Kalman filter requires 
knowledge of real time input forces for proper estimation. Therefore, the approach 
suggested in Figure 1-7 cannot be used to solve the rotor estimation problem. 

Another recently developed time domain Inverse System Identification method called the 
Sum of Weighted Accelerations Technique (SWAT) has been applied to many different 
impact problems [7]. The limitations of using SWAT lie in the fact that it can only 
reconstruct the sum of the external forces acting on a body’s center of mass and not the 
individual applied forces. To overcome this shortcoming, Genaro and Rade created a 
variation of SWAT which would yield the input forces [8]. However, this process 
introduces a shortcoming of its own in the fact that the number of sensors must be equal to 
or be greater than the number of responding modes, which is not the case for the CR. 

A further limitation of the input estimation or indirect force measurement techniques is that 
the process is found to be numerically ill-conditioned [9]. The numeric ill-conditioning 
occurs during calculations that require inverses of matrices which allows very small errors 
in measurements to result in large errors in estimated forces. 
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Another approach includes the use of general structured (GS) observers for state estimation 
during the case where inputs are unknown. A method for designing a full order unknown 
input observer (UIO) based on a GS observer is presented by Chang, You, and Hsu which 
allows for state estimation despite the existence of unknown inputs or uncertain 
disturbances [ 10]-[ 13]. However this method cannot be used for the CR problem since it 
requires the number of outputs (measurements) to be greater than the number of unknown 
inputs. For the CR problem, there are only 4 outputs versus the 8 possible inputs. 


1 .3 Thesis Overview and Content 

Chapter 2 provides a problem overview and a concise description of the CR system. This 
description includes a list of assumptions made during problem formulation and the process 
used to create a simplified model, which includes the rotor, the shroud, and a two-mass ISS 
flex model, for analysis and testing purposes. Also, the reference frames used as well as 
the derivations of the equations of motion for the simplified system are presented. 

Chapter 3 presents a detailed description of the proposed solution method by introducing a 
formulation of the disturbance models. The rotor disturbance is derived as a function of 
imbalance geometry, mass/inertia, and spin rate, while the ISS disturbance modeling is 
accomplished though a sinusoidal approximation of the effect of a pulse train through the 
system dynamics. Issues dealing with the peripheral effects of this sinusoidal 
approximation are examined along with both plant and filter system observability. 
Examples of both observable and unobservable filter models are presented using a modal 
form of the observability test. A time varying observability test is presented and used to 
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determine observability during cases of rotor spin-up and ISS operation. Also, the discrete 
Kalman filter equations and algorithm are introduced along with a method for calculating 
initial Kalman filter parameters. 

Chapter 4 provides a summary of the testing conducted to analyze estimation capabilities 
using the solution method proposed in Chapter 3. A description of the different 
performance measures used to evaluate Kalman filter performance is given. These 
measures included percent amplitude error in estimation, error covariance standard 
deviation envelope, error duration, and time to convergence. Testing was conducted to 
evaluate the validity of the proposed solution method and to show improved performance 
through disturbance model expansion within the filter dynamics. Finally, Monte Carlo 
analysis was performed to show both robustness of the estimator as well as its sensitivity to 
different uncertainties. The following computer programs were used to run all simulations 
and to perform all data analysis: Matlab Version 6.5.1.199709 (R13SP1) and Simulink 
Version 5.5.1 (R13SP1+). 

Chapter 5 provides a summary of the conclusions, along with a description of possible 
future work on the estimation process proposed in this thesis. 
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2 Problem Overview 

This section provides a problem overview and a concise description of the CR system. 
This description includes a list of assumptions made during problem formulation and the 
process used to create a simplified model, which includes the rotor, the shroud, and a two- 
mass ISS flex model, for analysis and testing purposes. Also, the reference frames used as 
well as the derivations of the equations of motion for the simplified system are presented. 


2.1 Modeling Assumptions 

For design and verifications purposes, a simplified model consisting of the CR system on a 
flexible ISS platform was used (see Figure 2-1). 


z 

A 


I 
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During the modeling process, in order to simplify the problem all external forces on the ISS 
other than attitude control jet firings, such as gravity gradients, ISS Control Momentum 
Gyroscopes (CMG) torques, aerodynamic forces, and orbital effects were neglected. 
Secondly, all masses which make up the ISS, shroud, and rotor are considered to be rigid 
bodies. The rotor is assumed to be cylindrical and therefore symmetric about the axis of 
rotation. ISS flexibility was modeled using a two mass-spring-damper system. A nominal 
ISS configuration was used for the determination of mass and inertia values of the two 
mass ISS flex model created for testing purposes. 

The origin of the reference frame is located at the geometric center, gc, which is the point 
on the x-y plane of the rotor through which all of the springs and dampers act. The gc is 
defined during equilibrium, and is the non-rotating inertial reference frame. 

The center of mass, cm, of the shroud and ISS flex model masses are collocated with the 
eg, thus eliminating any coupling between translation and rotation in or about the x or y- 
axis between the shroud and ISS flex model. Only the rotor’s static cm, noted on Figure 
2-1, is located directly above the reference frame along the z-axis, which causes coupling 
in the translational and rotational equations of motion between the rotor and the shroud. 


2.2 Derivation of the Linear, Time-Varying Equations of Motion 

In this section, the equations of motion for the simplified model has been developed to 
include the rotor, the shroud, and a two-mass ISS flex model. Each of the four masses has 
4 degrees of freedom (dofs) for a total of 16 dofs. 
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A cross-section of the model, in the x-z plane, that was used in the derivations of the 
equation of motion in x-axis is shown in Figure 2-2. This figure is not to scale. 



—JUT 


Rotational Spring 


Rotational Damper 


Figure 2-2. Simplified Model 


Typical ISS flex modes, between -0.01 and -1.0 Hz, are captured by a two-mass ISS flex 
model, which attaches to the shroud through translational and rotational springs and 
dampers. Each of the four masses in this simplified model has two translational dofs (x and 
y-axis) and two rotational dofs (about x and y-axis), resulting in a model with 16 dofs. In 
addition, the disturbance on the rotor, d ro u,r, caused by rodent motion during operation, acts 
on the rotor mass, while the disturbance on the ISS, diss, caused by jet firings, acts on the 
outside mass of the ISS flex model. Note that the relative measurement, x re i, is a relative 
measurement between the rotor ant the shroud. The system is time varying due to variable 
rotor spin rate experienced during rotor spin-up and spin-down. 
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2.2.1 EOMs for Coupled x Translation and <t> y Rotation 

The equations of motion describing the x translational motion and the <J) y rotational motion 
of the simplified model were calculated using Figure 2-3. 



K 1t> -^y Z^x 








Figure 2-3. Model Used for Derivation of X-translation and (jyrotation 


The rotational dofs shown in Figure 2-3 are relative to the inertial reference frame. 
Rotational and translational springs and dampers are located between each mass. Each 
mass is depicted separately to show translational and rotational dofs, but it is important to 
recognize that all masses of the same label are actually the same mass. That is to say that 
there is only one shroud (M s ), one ISS Mass 1 (Mi), and one ISS Mass 2 (M 2 ). This is also 
the case for the figure used for the derivation of y translational and (|> x rotational EOMs. 
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2.2.1 .1 EOMs for x Translation 


For the x translation equations of motion, 


F = nix 


( 2 - 1 ) 


where 


*r = 


M , 


x, = • 


M . 


x. = — — F. , , and 
1 M, u 


~^2 = 


M. 


2x 


( 2 - 2 ) 


The variable M represents each mass represented in the simplified model. The forces, F in 
equation ( 2-2 ), are provided by the springs, dampers, and external disturbance forces (See 
Figure 2-3). These forces are 


f rx = ~C Rtx x r + C Rlx x s + C^Lsin <p y - K Rlx x r + K Ru x s + K Rlx L sin <p y + d x 
F s< = C Ru x r -(C te +C Zu )x s +C Zlx x l -C Rlx Lsin <p y 

+ K Ru x r — (K Rlx + )x s + K Zlx x x — K Rtx Ls\n <p y ^ 

F \x = Cztx*s ~ (Cz« + C\ ix )*i + C Ux x 2 + F Zlx x s — ( K 2jr + K Ux )jr, + K Ux x 2 
F 2z = CfoX , — (C,„ + C 2tx )x 2 + K Ux x { — ( K Ux + K 2 u )x 2 + d x 


(2-3) 


where external disturbance forces are defined by Figure 2-4, 



Axis of Motion 
x = x translation 
Y = y translation 
<|> x = rotation about x 
<|> y = rotation about y 


K Location 
r = rotor 
s = ISS 


Figure 2-4. External Disturbance Naming Convention 


The naming conventions for the states, the stiffness/damping values, and the mass/inertia 
values can be found in Figure 2-5, Figure 2-6 and Figure 2-7, respectively. 
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J k 

* 

Axis of Motion ^ 

x = x translation 
Y = y translation 
0 X = rotation about x 
<j> y = rotation about y 


Which Mass 
r = rotor 
s = shroud 

1 = ISS mass 1 

2 = ISS mass 2 


Figure 2-5. State Naming Convention 


K = Spring 
C = Damper 

Location 
R = between Rotor and Shroud 
Z = between Shroud and ISS mass 1 

1 = between ISS mass 1 and ISS mass 2 

2 = between ISS mass 2 and inertial 



Direction of Motion 
x = in the x direction/about the 
y = in the y direction/about the 


Description of Motion 
t = translation 
r = rotation 


x axis 
y axis 


Figure 2-6. Naming Convention for Stiffness and Damping Values 



M = mass 
I = Inertia 


Which Mass ' or 
r = rotor 
s = shroud 

1 = ISS mass 1 

2 = ISS mass 2 


Which Inertia 
d r = rotor (transverse) 
p = rotor (axial) 
s = shroud 

1 = ISS mass 1 

2 = ISS mass 2 


Figure 2-7. Naming Convention for Mass and Inertia Values 
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Assuming small angle rotations allow for linear but coupled equation for x translation 
given in equation ( 2-4 ). Substituting equation ( 2-3 ) into equation ( 2-2 ) gives 


M r x r "*■ C r, x x r — C Rtx x s — C Rtx L tp V ' + K Rtx x r — K Ru x s — K Rlx L<p y = d Xr 
^ s — C Ru x r + (C Rlx + C /Jx )x s — C 7jx X\ + C Rtx L(jp y 

~ K Rtx x r ^ Rtx K- /Jx )x s — K Zjx x [ + K Rlx L<p >r = 0 

M x ‘x x — C Zlx x s + (C Zx + C Ux )x x — C Ux x 2 — K Zlx x s + ( K a;( + AT lw )jf| — K Xtx x 2 = 0 
M 2 x 2 — C Ux x | + (C Ux + C 2u )x 2 — K Ux x x + ( K Ux + K 2 ix )x 2 = d ^ 


( 2 - 4 ) 


2.2.1. 2 EOMs for <j> y Rotation 

The angular momentum equations, in the rotor body reference frame, provide the rotational 
equations of motion for the rotor are given by Yamamoto [15] [14] as 

T +r +/ „Wv, 

T«r=l d '<Py, 

7 * r = 0 

( 2 - 5 ) 


where l d is the transverse inertia for the rotor, / is the rotor spin axis inertia, a> is the 
rotor spin rate about the spin axis, and 

= Rty Ly r + C Rly Ly s —(C Rrx +C Rly L )<p x + C Rrx <p x ^ 

~ K R,y L y r + K R,y L y s ~ Rrx + K R,y L ~ )<Px r + K RrxPx, + d 0xr 

( 2 - 6 ) 


and 


r<?,r + C Rlx Lx r C Ru Lx s (C Rry + C Rlx L )<Py r + C Rry <Py t 


+ K Rlx Lx r K Rtx Lx s (K Rry + K Rix L )(p y + K RrvVv, 




( 2 - 7 ) 
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Since the rest of the masses are not involved in rotations, then H = Iw can be used along 
with the small angle approximation to calculate linear but coupled equation for and <t> y 
rotation 


A/A -lp(°<P Xr ~C Rlx Lx r +C Rlx Lx s +(C Rry +C Rtx L 2 )<p yr - C Rry ,<p y 

~ K Ru Lx r K Rix Lx s + ( K Rry + K Ru L )<p y — K Rry (p y> = d (jxr 

h<Py, " AA + <A + A A ~ A^, - A A +(^+^)^ ~^Zry<Py, =0 

An - An +(A + AA - An - An +<a + A A - An =0 
/ 2 n - An + <A + A A - An + (A + A>n = A 

( 2 - 8 ) 

The coupling occurs in the equation for the rotor motion due to the fact that the cm of the 
rotor is not collocated with the connection point of the springs/dampers. The distance, L, is 
used to account for translation in the x direction due to rotor rotation about the y axis at the 


cm and vice versa. 
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2.2.2 EOMs for Coupled y Translation and <|> x Rotation 


The EOMs describing the y translational motion and the <)> x rotational motion of the 
simplified model were calculated using Figure 2-8. 



Figure 2-8. Model Used for Derivation of Y-translation and (^-rotation 
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2.2.2.1 EOMs for y Translation 

For the y translation equations of motion, 


F = my 


where 


(2-9) 


^ = 7T F ’V’ y. =TT F ly' and =~Fr F 2y 

M s M , M 2 

( 2-10 ) 

M represents each mass of the simplified model. The forces, F, are provided by the 
springs, dampers, and external disturbance forces (See Figure 2-8). These forces are 



Pry = ~C R »y r + C Rty y s - C Rly L sin <p x - K Rly y r + K Rty y s - K Rty L sin (p x + d y> 
P sy ~ C Rlv y r — (C % + C ay )y s + C /Jv y t + C Rly L sin <p x * 

+ K R, y y r ~(^R,y + K Zty)ys + * 2**1 + K R<y L S™<P Xt 
P\y = ^ZryXv — (^Zry + y 1 + ^l/y^ + ^Zry^i — ^ Zty + P\ty )^l + ^llyV 2 
F 2 y = C lv y, - (C l0 , + C 2ty )y 2 + K Uy y, - (K Uy + K 2ty )y 2 + d y> 


(2-11) 


Assuming small angle rotations allow for linear but coupled equation for y translation 


M r y r + C Rly y r -C Rly y s + C Rty L<p x + K Rly y r - K Rly y s + K R , y L(p x = d y 
M s'y, -CR*y r + (C Rty +C Zty )y, -C^y, -C Rty L<p Xt 

— P Rn^r + (Priv + Pziy ) -Xs — P/Jy X I — P Rl^ j( P\, ~ 0 

+(C Z , V +C l(v )y l -C Xty y 2 -K 7ty y s +(^ + AT l(y )y, =° 

M 2^2 - C|, v y, + (C„ v + C 2fv )y 2 - K„ v y, + (K lty + A: 2/V )y 2 = d y 


( 2 - 12 ) 
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2. 2.2.2 EOMs for <J> X Rotation 

Since the rest of the masses are not involved in rotations, then H = Iw can be used along 
with the small angle approximation to calculate linear but coupled equation for and <|) x 
rotation 


+ +C Rly Ly r ~C Rly Ly s +(C Rrx + C Rly L 2 )<p x -C Rrx <p x 

+ K R„ L y r ~ K Rly Ly s + (K Rrx + K Rty L 2 )<p x - K Rrx <p X ' = d+ r 
7 A - C Rrx <P x , + (C Rrx + c zrx )<p x - - K Rrx <p x + (K Rrx + K Zrx )<p x - = 0 

I \P I, ~Czrx*Px, + ^Irjr _ ^Irx^x, ~^Zrx < Px, + ^ Zrx ^\rx )^r, ~^\rx < Px 2 = ^ 

h<Px 2 -C Ux <P Xl + (C,„ +C 2rx )<p Xi - K Wx <p Xi +(K lrx + K 2rx )<p X2 =d^ 

( 2 - 13 ) 


Again, the coupling between the y and <|> x occurs from the distance L between the rotor cm 
and the location where the springs and dampers attach. 

The time-varying aspect of the EOMs comes from the spin rate of the rotor, to, (seen in 
equations ( 2-8 ) and ( 2-13 )) which ramps up from 0 Hz to 0.7 Hz over a chosen time 
interval. Note that the z-axis translation of the VIM/CR is ignored since knowledge of the 
motion in the z direction is not necessary for control applications. 
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2.3 Model Frequency Response 

Values for the various model parameters such as mass/inertia and stiffness/damping values, 
were selected to capture the expected physical system dynamics. Using these assumed 
values, the frequency response from all 8 disturbance forces to the corresponding 4 relative 
measurements was calculated and is shown in Figure 2-9. 


Bode Plot: All Direct Transfer Functions 



Figure 2-9. Frequency Response of Reduced System 


The dashed lines show transfer functions including rotor disturbance inputs while the solid 
lines show transfer function including ISS disturbance inputs. 
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3 Solution Method 

During the estimation of the centrifuge rotor states, two main problems are encountered. 
First, the deterministic plant inputs (4 rotor and 4 ISS disturbance inputs) are not available 
to the observer for use in the estimation algorithm. Secondly, the ISS disturbance inputs 
into the plant are applied in the form of a pulse train of jet/thruster firings rather than as a 
sinusoid disturbance. This could pose a problem for the chosen solution method and will 
be discuss later in this Chapter. 

The first problem will be solved through the use of a plant model which has been 
augmented with disturbance states (See Sections 3.1.1 and 3.1.2). Estimation of the states 
for absolute rotor motion will be completed by using a Kalman Filter algorithm on this 
augmented plant model (See Section 3.4). The use of internal disturbance models will 
allow for estimation with a Kalman filter without the need for input measurements. This is 
vital to successful estimation, since normal Kalman filter operation requires knowledge of 
the inputs. 
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Not Viable -s 


Proposed 
Solution ^ 
Method 


Standard Kalman Filter Formulation 


d s 

d r 



Kalman Filtering with Internal Disturbance Model 

d. 


Plant 

X re , 

Observer 


+ 



Internal Disturbance Modeling 


A 

X 


Figure 3-1. Method Comparison with Standard Kalnian Filter Formulation 


Figure 3-1 shows the difference between the standard Kalman filter formulation, which 
explicitly includes the known, deterministic inputs, and the method that is employed in this 
thesis, which instead models these inputs as additional filter states, to circumvent the fact 
that the disturbances are not available as inputs into the observer. 

The second problem will be solved by using a sinusoidal approximation of the effect of the 
ISS pulse train disturbance on the plant dynamics, using only frequencies which result in a 
high gain through the plant. Frequencies where the plant attenuates the input signal are not 
important since they produce little effect on the measurement. Further explanation follows 
in Section 3. 1.2.1. 


27 


3.1 Filter Model EOMs 

The basic filter model, consisting of the rotor, VIM, and a 2-mass ISS flex model, has the 
same EOMs as the plant model (equations ( 2-4 ), ( 2-8 ), ( 2-12 ), and (2-13 )), but to 
allow for variation, the filter model coefficients will be allowed to deviate form the plant 
coefficient values. The filter model is signified by the addition of an T at the end of the 
coefficient variable names. In addition to modeling the plant within the filter, the rotor and 
ISS disturbances also need to be modeled. The process of integrating disturbance models 
into the filter model will be explained in Sections 3.1.1 and 3.1.2. 


3.1.1 Derivation of Rotor Imbalance Disturbance Forces 

The imbalance disturbance forces acting on the rotor have been derived as a function of the 
rodent mass, M rat , the transverse and axial rotor inertias, Ij r and I p , the distance of the center 
of mass (cm) from the spin axis, e, the spin rate, and the angle between the spin axis 
and the vector from the rotor tip to the rotor cm, a. 


$ 



(£>ri>£) " Rotating Coordinate Frame 

a - Phase angle from £ axis to the CM 
P - Angle from £ axis to the CM 
e - Distance from the origin to the CM 
in the plane 


Figure 3-2. Rotor-Fixed Rotational Reference Frame 
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The value for the scalar parameters e and a are used to defined a force vector in a rotor- 


fixed, rotating reference frame (r|, £,), then rotated via a coordinate transformation. 


equation (3-1 ), back into the inertial reference frame (x, y, z). 



cos# 

sin# 

V 

-sin# 

cos 6) 

c. 


( 3 - 1 ) 


where 0 equals G+t. The rotor imbalance disturbance equation in the inertial reference 
frame is listed below [15]. 


d xr = M£OJ 2 r Cos(w r t + p r ) d^ r = -(/, - l p )aw 2 Sin{(O r t + p r ) 
d yr = Meo))Sin{co r t + ) d^ r = ( / f - I p )aoj 2 Cos{a> r t + (3 r ) 

( 3 - 2 ) 

Note that these imbalance disturbance forces and torques found in equation ( 3-2 ), are not 
available to the filter as inputs. These disturbances are modeled as second order oscillators 
of the form 


^1 r Z 2r 


where the solutions to these differential equations are new states: z Xr and z 2r 

Z Xr = Cos(a) r t + fi r ) 
z lr = -co r Sin(co r t + p r ) 

( 3 - 4 ) 

Equations ( 3-2 ) can be rewritten in terms of these new disturbance states as 

d xr = Meo) 2 z lr d^ r =(I g -I p )a(o r z 2r 

dyr =-Mea) r z Xr d^ r = (/, -1 p )ao) 2 z lr 

( 3 - 5 ) 

Substituting these equations back into the equations of motion for the filter model results in 
a filter model with a state vector, jc, of length 34 and of the form 

x = [z zf z = [x r x s x x x 2 y r y s y x y 2 <j> x <f> x </> Xi <p Xi <t> y , <t> y , <t> h <t> yi z x J 
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3.1.2 ISS Disturbance Modeling 

The rotor imbalance disturbance is sinusoidal in nature due to rotor spin, while the ISS 
disturbance force is applied on the ISS in the form of jet impulses. The jet impulses result 
from the action of the ISS attitude control system [16]. The magnitude of the jet force is a 
constant, thus making the frequency of the jet firings and the overall on-time the only 
control variables. 


3.1. 2.1 Sinusoidal Approximation of a Pulse Train 

The effect of the ISS pulse train disturbance input on the output of the plant can be 
represented by a Fourier series and its related fundamental frequency. A sinusoid of that 
fundamental frequency can be used to model the effect of the pulse train in the observer 
model. An example is shown in Figure 3-3 where a pulse train input at 0.399 Hz is applied 
to the plant, and the output is a sinusoid with a frequency of 0.399 Hz. 


a) ds =0.399 Hz 


bW- 


Actual Jet Firing 
ISS Disturbances 



PLANT 


Vel 


• Output From Pulse Input 


PSD of Output from Pulse Input 


A 




.MO Plot of K ff| Output duo to Mm Input 



Figure 3-3. PSD of Plant Output due to Pulse Train Input is Sinusoidal 
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If the ISS pulse train disturbance excites the plant at a high gain frequency, then the output 
will have a large contribution as a result of the ISS disturbance. However, if the ISS pulse 
train disturbance excites the plant at a low gain frequency, then the output will have a small 
contribution resulting from the ISS disturbance. Therefore, it is important to determine the 
plant peak gain frequencies. These frequencies are determined from the frequency 
response plots (see Figure 3-4). These peak frequencies will now be used for ISS 
disturbance modeling within the filter. The discussion of the effect of the amplitude 
mismatch between the sinusoid and the jet firing will be conducted in Section 3. 1.2.4. 


Bode Plot: All Direct Transfer Functions 



Figure 3-4. 


Frequency (Hz) 

Frequency Response Plots 
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3.1. 2.2 ISS Model 1: 4 ISS Disturbance States 

Using the approach discussed in Section 3. 1.2.1, an ISS disturbance model can be created 

representing the jet firing as a sinusoidal disturbance with single translational and single 

rotational disturbance frequency of CDdst and 0)d Sr , respectively, of the form: 

d xs = F,Cos(co dJ + 0 S ) d 0s = -T*Sin(a) d j + y s ) 
d ys = F y Sin(ci) ds ,t + 0 S ) d^ = T^Cos(co d j + y s ) 

( 3 - 6 ) 

The ISS disturbance amplitudes F x , F y , T^, and T^y, are the known specifications of the 
thrusters located on the ISS. The ISS disturbance forces found in equations ( 3-6 ) are 
modeled as second order oscillators of the form 

^1j = ^3i 
^2s ^4j 

^3 i = Z 1 j 

^4.v = ~ a) dsr Z-2s 

( 3 - 7 ) 

where the solutions to this set of differential equations are the new states: z l5 , z 2s , z, s , and 
Z 2s , where 

Zi. = Cos(a*,t + B s ) z u = -co dsl Sin(co d j + B s ) 

z 2s = Cos(a) d j + y s ) z 2s = -co dsr Sin((o dsr t + y s ) 

( 3 - 8 ) 

Using equations ( 3-8 ) the ISS disturbance equations can be rewritten in term of the new 
ISS disturbance states as 


d xs F X Z\ S d 4ks 


0X • 


CO, 


t Isr 


. . _ T 

“ ys Ms I «y^2s 


CO, 


,kl 


( 3 - 9 ) 

Equations ( 3-9 ) can be substituted back into the equations of motion for the filter model. 
The resulting filter model has a state vector, jc, of length 38 of the form 
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x = [z if z = [x r x s jc, * 2 y r y s y, y 2 <p Xr <j> x ^ </> Xi ^ <z> y , <p h z u z lr z u z 2s ] 

This ISS disturbance model will be called “ISS Model 1” for testing purposes. 


3.1. 2.3 ISS Model 2: 8 ISS Disturbance States 

ISS Model 1 is expanded to create a new ISS disturbance model labeled “ISS Model 2”, 
which contains 8 ISS disturbance states for a total of 12 disturbance state, when including 
the rotor disturbance states. ISS Model 2 represents the ISS disturbance as sums of two 
sinusoidal disturbances, with two translational and two rotational disturbance frequencies 
of (Odsti. Wdsi 2 , and (iWi, GW 2 , respectively, of the form: 

d »\ = F xiC° s ( a >d S ,i t + A)+ F x iCos(0)^ 2 t + P s ) = -T^Sinico^t + V x )-T« 2 Sin(co <lsr2 t + y s ) 

d y *\ = F yi Sin(co dul t + p s ) + F y2 Sin(eo (ls2 lt + p s ) d* sl = T^Cos(co dirl t + y s ) + T^ 2 Cos(co ihr2 t + y s ) 

( 3 - 10 ) 


The two ISS disturbance amplitudes in each axis are assumed to equal the known 
specifications of the thrusters on the ISS (i.e. F x = F X | = F X 2 , F y = F y i = F y 2 , T$ x = T^i = 
T(j) X 2 , and T$ y = T$ y i = T^). The ISS disturbance forces found in equations ( 3-10 ) are 
modeled as second order oscillators of the form 

Z|j = Z 5 j Z 5j = ~ (t) dst\ Z 1j 

?2 s = Z 6j Z 6 s = ~^dsl2 %2s 

%3s ~ Z 7j ^7 s ~ ~ ( °dsr\ ^3i 

^4j = Z 8j Zg s = ~ (d dsr2 Z*s 

( 3 - 11 ) 

where the solutions to this set of differential equations are the new states: , z 2s , z 3s , 


> z Ix > z 2s • z 3j - and z 4s where 
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*i, = CosiQi^t + # ) z u = -(O dsll Sin(a) dsn t + fi t ) 
z 2 s = Cosito^t + A ) z 2j . = -(0 dsl2 Sin((0 dsl2 t + fi s ) 
z 3s = Cos(co Asrl t + y s ) z 3s = -a) dsrl Sin((O dsrl t + y s ) 
z a.s = Cos(o) dsr2 t + y s ) z 4s = ~(O dsr2 Sin((O dsr2 t + y s ) 

(3-12) 

The ISS force and torque inputs ( 3-10 ) can be rewritten in terms of the ISS disturbance 
states as 


d xs = F x\Zu + F x2Z 2s 


0™ 


T T 

. t 1 <fix2 . 

^3 s m *4, 


d — ^ 


dsrl 


d =-— z, -- F>2 

}8 ^ ^lS 


at*. 


dsl 1 




^2s d <9v.\ F ^y\Z'is ^2 *4, 


*1 2 


(3-13) 

Equations ( 3-13 ) can be substituted back into the equations of motion for the observer 
model, increasing the number of states from 34 to 42. The resulting filter has a state 
vector, x, of length 42 and is of the form 



z = k jc, •«. * 2 y r y s y, y 2 K <t>*. K K K K 0 y , K 


Z\ r r ^ls ^2 s ^3j ^4 


3.1. 2.4 PSD Difference between Sinusoid and Pulse Train 

It is important to note that there is a power spectral density (PSD) difference in plant 
measurements (outputs) between a sinusoid disturbance and a pulse train disturbance; a 
sinusoidal disturbance of the same frequency as the pulse train with a small on-time creates 
an output with a higher PSD. The Kalman Filter determines the amplitude of the modeled 
disturbance sinusoid in order to get an equivalent sinusoid which would have created the 
same output as that from the pulse train input. Therefore, it is necessary to take into 
consideration the difference in the PSD during ISS disturbance state comparisons. Figure 
3-5 illustrates the difference in the PSD of plant outputs when excited by a sinusoidal force 
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and a pulse train input force. The pulse on-time, defined as the duration per cycle for 
which the value is not equal to zero, is 20% of the period, and the frequencies of the 
sinusoid and the pulse train are the same. 


Modeled Sinusoidal 
ISS Disturbances 

k/V¥V-> 


O) ds =0.399 Hz 
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Figure 3-5. Output PSD Differential Between and Sine and Pulse Train Inputs 


The power spectrum is generated by using a Fourier transform and taking the square of the 
magnitudes of the complex coefficients [17]. Therefore, the effect of using a pulse train 
rather than a sign wave can be calculated as: 


PSD Ratio = 



V^sine 


(3-14) 

where P S j ne is the PSD due to the sine input at the excitation frequency (0.399 Hz in the 
example) and P pu i S e is the PSD due to the pulse train input at the same excitation frequency. 
For the example presented in Figure 3-5, the PSD ratio equals: 


PSD Ratio = 



V 2.6435*10 -J =ft3935 
VI .7073* 10 1 


This result can be interpreted as the factor by which the original sinusoidal input amplitude 
would need to be multiplied by in order to get the equivalent plant output with a pulse train 
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input. By multiplying the original sinusoid by 0.3935, it can be shown that the PSD of the 
output is now exactly the same. 


_ • _ PSD of Output from Sine Input 



Figure 3-6. Output PSD Results of Amplitude Ratio Sine and Pulse Train Inputs 


Since the ISS disturbance frequency could be any value between 0.01 Hz to 1 Hz, the 
entire range of frequencies was scanned to determine the actual equivalent disturbance PSD 
ratio for all ISS disturbance inputs. 


PSD Ratio (x) 


Average 



Fraquancy (Hz) 

Figure 3-7. PSD Ratio Over All ISS Frequencies for All ISS Disturbance Inputs 
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The average PSD ratio over the range of possible ISS frequencies is equal to -0.3933 for 
each disturbance input. Therefore, when comparing the disturbance state estimate to the 
actual disturbance state, the PSD ratio will have to be factored into the actual disturbance 
state amplitude. This is helpful in determining estimation performance. 


3.2 Observability 

Observability of the system is necessary in determining the viability of the Kalman filter as 
a solution method. The available methods used for determining plant observability include 
the well known Popov-Belevith-Hautus (PBH) Criterion [18]-[20] as well as a modal 
criterion for observability [21]. Both methods will be examined, but due to problems with 
ill-conditioning, the modal criterion for observability will be used to determine 
observability of the plant and filter models under time invariant conditions. A time varying 
observability test will be used to determine observability during operations such as rotor 
spin-up or spin-down and ISS maneuvering. 


37 


3.2.1 PBH Criterion for Observability 

Consider a continuous time system described by 


x = Ax + Bu 

( 3-15 ) 


y = Cx + Du 


( 3 - 16 ) 


where x = state vector (n - vector) 

y = output vector (m - vector) 

A = System Dynamics (n x n matrix) 

B = Input Matrix (n x r matrix) 

C = Output Matrix (m x n matrix) 

D = Direct Transmission Matrix (m x r matrix) 


The solution to equation (3-15) is 


x(t) = e Al x(0)+ je A<, ' r, Bu(r)dr 
0 


( 3 - 17 ) 


and y(t) is 

t 

y(t) = Ce Al x(0) + C je A<t ' r) Bu(r)d r + Du 

0 

( 3 - 18 ) 

Since the matrices A, B, C, and D are known and u(t) is also known, then the last two terms 
on the right half side of equation ( 3-18 ) are known quantities. Therefore, they can be 
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subtracted from the observed value of y(t), and it is sufficient to consider the unforced 
system described by 


x = Ax 


( 3 - 19 ) 


y = Cx 

Referring back to equations ( 3-17 ) and ( 3-18 ) 


( 3 - 20 ) 


x(t) = e A, x(0) = 2> k (t)A k x(0) 

k=0 


( 3 - 21 ) 


and y(t) is 


n-1 

y(t) = Ce Al x(0) = £tf k 0)C A k x (°> 


k=0 


( 3 - 22 ) 


For the system to be observable, given the output y(t) over a time interval 0 < t < t|, x(0) is 
uniquely determined from equation ( 3-22 ). It has been shown that for this to occur, the 
rank of the Observability matrix, O, of size (n x nm) must be full (i.e rank(O) = n). This is 
the so called PBH criterion for observability [18][19], 


0 = [c CA ••• CA"- | J r 

( 3 - 23 ) 

The problems with using the PBH criterion occur if some eigenvalues of A are greater than 
one while others are less than 1. Since the observability matrix, equation ( 3-23 ), requires 
A"' 1 , then if the number of state, n, is large, then the observability matrix will become 
numerically ill-conditioned. The singular values less than 1 will trend towards zero while 
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the singular values greater than 1 become large and the value of n increases. Since the rank 
test is determined by the number of singular values above a certain tolerance (10' K> for 
Matlab rank command), as n approaches 32, the number of singular values which fall 
below the tolerance increases. Since the observability matrix requires A n_l (A 31 when using 
the plant dynamics), the rank of the observability matrix, is only 11; a rank of 32 is 
required for full rank. 

The condition number, which is used to measure the level of ill-conditioning, is defined as 
the ratio of the maximum singular value to the minimum singular value. The larger the 
condition number, the more ill-conditioned the problem becomes. The observability matrix 
has a condition number of 3.432 x 10 24 . This shows that with the parameters chosen for 
the example, a severe problem of ill-conditioning does exists. Therefore, an alternate 
method is needed to determine observability. The modal criterion for observability 
eliminates the need to compute high powers of the system dynamics. 


3.2.2 Modal Criterion for Observability 


The modal criterion for observability is described by Ogata [21]. Consider a system 
described by equations ( 3-19 ) and ( 3-20 ). Also suppose that the A matrix is 
diagonalizable with the use of a transformation matrix, T, such that 


T"'AT = A 


where A is a diagonal matrix. Let us define 


( 3 - 24 ) 


x =Tz 


( 3 - 25 ) 
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where z is the transformed state. In terms of the new transformed states, equations ( 3-19 ) 
and ( 3-20 ) become 


z = T 'ATz = Az 


( 3 - 26 ) 


y = CTz 


( 3 - 27 ) 


Using equations ( 3-21 ) and ( 3-22 ) 

y(t) = CTe*‘z(0) 


( 3 - 28 ) 


or 


y t <«) = S(clT i )e 1 -'z 1 (0) 

i=l 

( 3 - 29 ) 

where n equals the number of states, c,[ denotes the k lh row in the C matrix, and X\ denotes 
the i lh eigenvalue. If cjTj =0 then the i th mode is unobservable in the k ,h output. If 
CT| = 0 then the i ,h mode is unobservable from all outputs. 


In other words, the system is observable if none of the columns of the m x n matrix CT 
consist of all zero elements. This is easy to see, since with the decoupled dynamics, if the 
i th column is found to be all zeros, then the corresponding state Zj(0) will not be a part of the 
output equation. 

If the system includes complex conjugate eigenvalues, then a modal A matrix can be 
created where the real eigenvalues appear on the diagonal of the matrix and the complex 
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conjugate eigenvalues appear in 2-by-2 blocks on the diagonal of the modal A matrix. For 
example, a system with eigenvalues (A.|, h. A*), the modal A matrix is of the form 



A 

0 

0 

0 


0 

a 

(O 

0 

A = 

0 

(O 

<7 

0 


0 

0 

0 

a 2 


( 3-30 ) 


where a = Re(A. 2 ) and co = Im(A. 2 ). 


To test observability for complex conjugate eigenvalues, both the real and imaginary dot 
products must be zero for that mode to be unobservable. In this case, columns of all zeros 
in the CT matrix would come as single columns for real eigenvalues and as adjoining 
columns for complex conjugate pairs. 

The modal condition for complete observability is also useful because using the inverse of 
the transform makes it possible to determine the combination of original states which cause 
the transformed state to be unobservable. 

z = T'x 

( 3-31 ) 

After determining which z states are unobservable by using the CT matrix, it is possible to 
determine the combination of x states which make up those unobservable z states by 
looking at the corresponding rows in the T 1 matrix. 
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3.2.3 Time Invariant Observability Test for Time Varying Dynamics 

It is important to check observability of both the plant dynamics and the filter dynamics. If 
the plant is not observable, then the filter will be unable to produce accurate estimates of 
the plant states, even if the filter model is observable. On the other hand, if the filter model 
is not observable, then even if the plant is observable, the filter will be unable to produce 
accurate estimates. Therefore it is important to check the observability of both plant and 
filter models. Testing could require any combination of the time varying parameters (rotor 
spin rate and ISS disturbance frequencies) to be held as a constant value. Therefore, some 
assurance of observability is required over all possible combinations of time varying 
parameters. 


3.2.3.1 Observability of Plant Model 

The time varying components of the plant dynamics are simply functions of the rotor spin 
rate. In order to test the observability for all spin rates the modal observability test was 
performed on the plant dynamics with spin rates from 0.001 to 1 Hz, with a frequency step 
size of 0.001 Hz. If any of the columns of the CT matrix are all zeros (or less than a 
tolerance of 10 5 ) then the system is considered unobservable. A observability plot (Figure 
3-8), where 1 means that the system is observable at the given spin frequency, shows that 
the plant is observable for the entire range of spin rates. 
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Observability Test (over all j 



Figure 3-8. Plant Observability Test over the Entire Range of Spin Frequencies 


3.2.3.2 Observability of Filter Model 

The time varying components of the filter dynamics are simply functions of the rotor spin 
rate and the ISS disturbance frequencies. ISS disturbance frequencies were assumed to 
have a range of 0 to 1.2 Hz and were discritized with a frequency step size of 0.001 Hz. In 
order to test the observability of the filter model, the modal observability test was 
computed using the filter dynamics for every possible combination of time varying terms 
(i.e. each disturbance frequency was tested for a given spin frequency). If any of the 
columns of the CT matrix were all zeros (or less than a tolerance of 10' 5 ) then the system 
would be considered unobservable. Figure 3-9 below shows that filter dynamics are 
unobservable only during extremely low spin frequencies (below 0.015 Hz). This is not of 
concern for two main reasons: 1) At such a low spin frequency, the rotor imbalance force 
will be negligible and 2) Since the rotor spins up from 0 to 0.7 Hz, the spin frequency will 
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be under 0.015 Hz for less than 7 seconds, assuming a 300 second ramp up period, 
therefore the system will only be unobservable for a very short duration (see Figure 3-9). 


Observability Test Results (Binary) 



Rotor Spin Frequency (Hz) 

Figure 3-9. Observability Test for All Time-Varying Dynamics Combinations 


3.2.3.3 Example of Fully Observable Filter Model 

The following are the results of the modal observability test for a fully observable case (to r 
= 0.7 Hz and all Q)d S = 0.4 Hz). The 1-norm of each column of the CT matrix was taken to 
simplify data interpretation. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Eigenvalues of 
A matrix 

-0 11682 + 
10 424i 

-0.11682- 
10 4241 

-0 073389 + 
1 0.31 7i 

-0 073389 - 
1 0 31 7i 

-0.18089 + 
7.6069i 

-0 18089- 
7 6069i 

-7 3227e-5 + 
6 5262i 

-7 3227e-5- 
6 5262) 

-2 3473e-5 + 
4 9663) 

-2 3473e-5 - 
4 9663i 

Norm of 
Columns of CT 
(w/Tolerance) 

0.09188 

0.09210 

0.07381 

0.07412 

0.08462 

0.08461 

0.00184 

0.01591 

0.01150 

0.02006 






11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

Eigenvalues of 
A matrix 

-0 0014442 + 
4 3275 

-0 0014442 - 
4 3275 

-0 062997 + 
3.5895 

-0 062997 - 
3 5895i 

-0 063107 + 
3.5478) 

-0 063107 - 
3 5478i 

-0 0026649 + 
2 5059i 

-0 0026649 - 
2 5059i 

-1 877e-5 + 
3 7725 

-1 877e-5 - 
3 7725 

Norm of 
Columns of CT 
(w/T olerance) 

0.01318 

0.19392 

0.03083 

0.31019 

0.04175 

0.29286 

0.04339 

0.37691 

0.00210 

0.00364 











21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

Eigenvalues of 
A matrix 

-0 043358 + 
1 .61 9i 

-0 043358 - 
1 61 9i 

-0 0046068 + 
1 391 3i 

-0 0046068 - 
1 3913i 

-000016548 + 
0 38928) 

-0 00016548- 
0 38928) 

-0 036858 + 
078154) 

-0 036858 - 
0 78154) 

-0 033126 + 
0 76038) 

-0033126- 
0 76038) 

Norm of 
Columns of CT 
(w/Tolerance) 

0.40752 

0.40720 

0.42679 

0.37970 

0.44470 

0.19488 

0.71313 

0.37802 

0.76313 

0.41044 







31 

32 

33 

34 

35 

36 

37 

38 


Eigenvalues of 
A matrix 

-0 0012268 + 
0 63498i 

-0 0012268 - 
0 634981 

5.551 1e-17 + 
3 7699i 

5 5511e-17 - 
3 7699i 

7.0777e-16 + 
0 6346i 

7 0777e-16 - 
0.6346i 

0 + 4 3982) 

0 - 4 3982) 

Norm of 
Columns of CT 
(w/Tolerance) 

0.50741 

0.07915 

0.00007 

0.00003 

0.37055 

0.05411 

0.00119 

0.00118 


Table 3-1. 1- Norm of the Columns of the Modal Observability Matrix (CT) 
for a Fully Observable Case 


Any number in the CT matrix that is less than a tolerance value of 10 5 was set to zero. The 
1-norm of the rows of the CT matrix, shown in Table 3-1, contains no zero values and is 
therefore fully observable. 

3.2.3.4 Example of Unobservable Filter Model 

The following are the results of the modal observability test for an unobservable case (o) r = 
0.01 Hz and all 0)ds = 0.7 Hz). Similar to the observable test case, the 1-norm of each 
column of the CT matrix was taken to simplify data interpretation. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Eigenvalues of A 
matrix 

-0 0805 ♦ 
10 336i 

-0 0805- 
10 338i 

-0 0801 ♦ 
10.3351 

-0.0801 - 
10 3351 

-7 38O-005 ♦ 
6 52621 

-7 38O-005 - 
6 5262i 

-3 33O-005 4 
4 96631 

-3 33O-005 - 
4.96631 

-0 00149 4 
432781 

-0 00149- 
4 32781 

Norm of Columns 
of CT 

(w/T olerance) 

0.06352 

0.08756 

0.06335 

0.08746 

0.00076 

0.01465 

0.00364 

0.04355 

0.01654 

0.18903 



11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

Eigenvalues of A 
matrix 

-0 00213 ♦ 
2 50511 

-0 00213- 
2 50511 

-0 0977 4 
3.37921 

-0.0977- 
3 37921 

-0 0663 4 
3 44231 

-0 0663 - 
3 44231 

-0 128 4 3 64291 

-0 128 - 3 84291 

-0 0906 4 
3 62351 

-0 0906 - 
3 62351 

Norm of Columns 
of CT 

(w/Tolerance) 

0.01778 

0.35617 

0.13410 

0.22749 

0.10347 

0.24106 

0.10366 

0.24204 

0.05677 

0.26287 



21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

Eigenvalues of A 
matrix 

-0 000217 ♦ 
1 41121 

-0 000217- 
1 41121 

-2 07e -005 ♦ 
3.77251 

-2 07O-005 - 
3.7725i 

-0.0351 4 
0 76881 

-0 0351 - 
0 76861 

-0 0364 4 
0 780161 

-0 0364- 
0 780161 

-0 00014 4 
0.390941 

-0 00014 • 
0 39094t 

Norm of Columns 
of CT 

(w n olerance) 

0.01993 

0.07976 

0.01219 

0.01677 

0.73708 

0.00972 

0.75301 

0.01150 

0.40435 

0.00686 












31 

32 

33 

34 

35 

36 

37 

38 

Eigenvalues of A 
matrix 

-0 0011 ♦ 
0 635061 

-0.001 1 • 
0 635061 

6 24e-017 4 
3.76991 

6 240-017- 
3 76991 

2 440-016 4 
0.0628321 

2 440-016- 
0.0628321 

1 380-016 4 
0.63461 

1 380-016- 
0 63461 

Norm of Columns 
of CT 

(w/T olerance) 

0.45802 

0.03761 

0.00032 

0.00024 

0 

0 

0.33446 

0.04833 

Table 3-2. Mod 
(CT) 

al 1- Norm of the Columns o; 
i for an Unobservable Case 

the Modal Observability Matrix 


Also, the same 10 5 tolerance was used in order to determine zero values. The table 
including the 1-norm of the columns of the CT matrix, shown in Table 3-2, does contain 
values of zero. Therefore, zero column vectors exist in the CT matrix, resulting in 
transformed states 37 through 40 being unobservable. Equation ( 3-31 ) is used to solve for 
the x state combinations which result in the unobservable transformed states. The rows of 
the T' 1 matrix which account for the unobservable transformed states can be found in Table 
3-3 below. The x state combinations which make up the z states can be calculated using 


Table 3-3. 
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x state 

x states 
1-16 

17 

x states 
18-35 

36 

x states 
37,38 

Zlr 

Zlr 

z state 

Z 35 

0 

1.002 

0 

0 

0 

Z 36 

0 

15.947 


Table 3 


-3. Rows of the T' 1 Matrix which Account For Unobservable Z States 


Therefore, 


z 35 — 1.002z, r 
z 36 =15.947z lr 


This shows that none of the rotor disturbance states can be observed at low spin rates. This 
finding is reasonable because since the rotor disturbance amplitude is a function of the 
square of the spin rate,<y 2 , a small spin rate would equal a very small disturbance 
amplitude; essentially, there is nothing to be observed. 
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3.2.4 Time Variant Observability Test 

Under conditions of spin-up or ISS maneuvering, the time invariant observability test is not 
sufficient in determining observability. An observability test is required which takes into 
account the time-varying system dynamics. Gelb discusses observability under the 
assumption of a time-invariant system [23], however, his approach applies to time-varying 
systems as well. Consider the discrete system 


x = Ax, x(t 0 ) = x 0 

0 - 32 ) 

II 

o 

X 

( 3 - 33 ) 

The solution to ( 3-32 ) is 

x(t) = <l>(t,t 0 )x 0 

( 3 - 34 ) 

where O(t,t 0 ) is the solution to the matrix differential equation 

-j-(4>(t,t 0 ))= A(/)4>(t,t 0 ) 
at 

( 3 - 35 ) 

where 

®(t 0 . ^0 ) = I 

( 3 - 36 ) 
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Using equations ( 3-33 ) and ( 3-35 ) 

yfo ) = c ('o ) <1 >(t 0 . t 0 )x 0 = C(t 0 )x 0 

^.) = C(t,)<I>(t 1 ,t 0 )x 0 

y(/ 2 ) = C(r 2 )<D(t 2 ,t 0 )x 0 

y(fn - 1 ) = )^ > (t n -l ’^0 ) X 0 

(3-37) 


’ y(0" 


c(t,) 

y(t.) 

= 

c(t,)o(t„t 0 ) 

_y(t 0 -i)_ 


C(t 2 )O(t n ,,t 0 )_ 


Z:n*m x n 

( 3-38 ) 


where n is the number of states and m is the number of measurements. The condition for 
which xo is observable for the measurement times to, ti, t n -i is that 


rank(Z) = n 


(3-39) 


Matlab and Simulink can be used to numerically calculate O(t,t 0 ), using the following 


algorithm 



> ( K'o’'o) 

'o) 


Figure 3-10. Numerical Solution for O(t,to) 
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where A(t) is the time varying dynamics and the <t> matrices are captured at each time step, 
i. Using the algorithm described in Figure 3-10 and equation ( 3-38 ), the plant, the filter 
using ISS Model 1, and the filter using ISS Model 2 are all full rank and therefore 
conditionally observable. 


3.3 Introduction to Optimal Linear Filtering 

Now that it has been proven that the system is observable, a filter can be used to estimate 
the necessary rotor states. The term filter refers to the estimation of state at the present 
moment using previous measurements. An unbiased estimate is one whose expected value 
is the same as the expected value of the quantity being estimated. A minimum variance 
estimate has the property that its error variance is less than or equal to that of any other 
unbiased estimate. A consistent estimate is one which converges to the true value as the 
number of measurement increase. By these definitions provided by Gelb [23], we are 
looking for an unbiased, minimum variance, consistent filter. 

When a controller requires state feedback, but the available measurements do not include 
all necessary states, there must be a method of estimating the missing states that contains 
minimal error. This requires the following [24]: 

• The ability to define a state-estimate error metric to be minimized in estimation 

• A knowledge of measurement error statistics, dynamic system models, and system 
input statistics 

• Algorithms for using this information to compute minimum-error state estimates 
The Kalman Filter is one such suitable algorithm for state estimation. 
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3.4 Advanced Optimal Linear Filtering: The Kalman Filter 

In 1960, R. E. Kalman published his paper entitled, “A New Approach to Linear Filtering 
and Prediction Problems,” describing the use of a recursive filter to solve the Weiner 
problem for gauss-markov sequences through the use of state-space representation from the 
viewpoint of conditional distributions and expectations [25]. The Kalman filter is powerful 
because it not only supports estimation of the past, present, and future, but can do so even 
if the modeled system is not known precisely. The following concise derivation is from 
Welch and Bishop [26]. 


A discrete time process, with a state vector xg91" is governed by the following linear 
stochastic difference equation 


x k = A k .,x k _, +Bu k _, + w k _, 

with a measurement output vector ze ST" 

z k =Cx k +v k 


( 3 - 40 ) 


( 3-41 ) 


where w k and v k are process noise and measurement noise, respectively. These noise terms 
are assumed to be white noise, having a Gaussian distribution with a mean of zero and a 
covariance of Q and R, where Q is the process noise covariance and R is the measurement 
noise covariance. 

p(w)~ N(0,Q) 
p(v)~ N(0,R) 

( 3 - 42 ) 
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The A, B, and C matrices from equations ( 3-40 ) and ( 3-41 ) are standard discrete-time 
state space matrices, with the following dimensions: Ae9V lxn , Be9^ nxr , Ce 9l mxn . The 
input vector has the dimension u€ 5H r . 

Given that 


x' k = a priori state estimate at time step k, given knowledge of the 
process prior to step k 

x k = a posteriori state estimate at step k given measurement zy 
where the a priori, ek\ and a posteriori, e*, state estimate errors are 


e 


k 


= X, 


-x v 


e k 


= X, 


-X 


k 


( 3 - 43 ) 


Therefore the a priori estimate error covariance, Pk", and the a posteriori estimate error 
covariance, Pk, are 

P k = E[e k e k ] P k = E[e k e k ] 


( 3 - 44 ) 


Also, the Kalman filter works in such a way that it sets the estimated state’s value at the 
expected value of the actual state. 


x k =E[x k ] 


( 3 - 45 ) 


Therefore, it is important to point out that the Kalman filter maintains the first two 
moments of the state distribution found in equation ( 3-44 ) and equation ( 3-45 ) above. 
Because of this, the a posteriori state estimate, x k , reflects the mean of the state distribution 
and the a posteriori estimate error covariance reflects the variance of the state distribution 
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[23]. In other words, the conditional probability density function of x k , conditioned on a 
value of z k , is defined as 

P(x k |z k )~ N(E[x k ],E[(x k -x k )(x k -x k )' ]) 

= N(x k ,P k ) 


( 3 - 46 ) 


To derive the Kalman filter equations we begin with the goal of finding an equation that 

A A _ 

computes X k as a linear combination of x k and a weighted difference between the actual 


measurement, Z k , and a measurement prediction Cx k . 


* k =x k +K k (z k -Cx k ) 

( 3 - 47 ) 

where the term multiplied by the gain, K k , is the measurement residual. If the measurement 
residual is equal to zero (the difference between estimated output and actual output is zero) 
then the a priori state estimate does not need to be altered before it becomes the state 
estimate at time step k. The n x m matrix K k , called the Kalman Gain, is used to minimize 
the a posteriori error covariance, P k , found in equation ( 3-44 ). 


Given equation ( 3-43 ) and by using equation ( 3-47 ) you get 



Substituting equation ( 3-48 ) into equation ( 3-44 ) you get 


( 3 - 48 ) 


P k = E[(x k - x k -K k (z k -C k x k ))(x k -x k -K k (z k -C.x,)) 1 ] 


( 3 - 49 ) 
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After performing the indicated expectation, then taking the derivative of the sum of the 
diagonal terms of the result with respect to k, and setting that result equal to zero, the 
equation can be solved for the Kalman Gain. 


K. =P l C;(c i P k C[ + R)‘ l 


( 3 - 50 ) 


A more rigorous derivation of the Kalman gain is provided by Gelb [23] and by Mangoubi 
[22] and follows below. Gelb begins with the following assumed form of a linear, 
recursive estimator 


= K k x k +K k z k 


( 3 - 51 ) 


where K k and K k are time-varying weighting matrices to be defined later. 

Given that 

x k =x k +e k 
x k =x k +e k 

( 3 - 52 ) 

by substituting equations ( 3-41 ) and ( 3-52 ) into equation ( 3-51 ) results in the following 
definition for a posteriori error at time step k. 

e k = (^k "*"K k C k -I]x k + K k e k + K k v k 

( 3 - 53 ) 

Since v k is defined as white noise with a mean of zero, the expectation of v k = 0. Then, if 
the expectation of the a priori estimation error equals zero (E[ek] = 0) then the estimator is 
unbiased (i.e. E[e k ] = 0) for any state vector if the bracketed terms in equation ( 3-53 ) are 
equal to zero. Therefore 
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[K k +K k C k -I] = 0 


( 3-54 ) 


and K k must be defined as 


K k =I-K k C k 


( 3-55 ) 


Rearranging the terms in equation ( 3-52 ) and substituting them into equation ( 3-53 ) 
results in 


*k -x k — [I"K k C k + K k C k + I]x k +[I K k C k ](x k x k ) + K k v k 
x k =[I-K k C k ](x k -x k ) + K k v k +x k 

= [I — K k C k ]x k — [I — K k C k ]x k + K k v k +x k 
= [I-K k C k ]x k -x k + K k C k x k + K k v k + x k 
= [I — K k C k ]x k + K k [C k x k + v k ] 

( 3-56 ) 


By substituting equation ( 3-41 ) into equation ( 3-56 ) the state update is obtained. 


x k =[I-K k C k ]x k +K k z k 

or 


Xk — x k + K k [z k -C k x k ] 


(3-57) 


Using equation ( 3-41 ), ( 3-52 ), and ( 3-57 ) the error dynamics are 

e k = ■ Kk^k ]e k +K k v k 


( 3-58 ) 


This equation for ek is used in order to update the error covariance Pk defined in equation 
(3-44). 
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P k = E[((I - K k C k )e k + K k v k X(I - K k C k )e k + K k v k )] 


( 3 - 59 ) 


Expansion of this equation leads to 

P k = E[(I - K k C k )e k e k ' (I - K k C k ) T + (I - K k C k )e k v* K k 
+ K k v k e k T (I-K k C k ) T +K k v k v k K k ] 


( 3 - 60 ) 


By definition 


and 


E[e k e k T ] = P k 


( 3 - 61 ) 


E[v k v[] = R k 


( 3 - 62 ) 


Since measurement errors are uncorrelated 

E[e k v[] = E[v k e k T ] = 0 

( 3-63 ) 

By substituting equations ( 3-61 ), ( 3-62 ), and ( 3-63 ) into equation ( 3-60 ) the error 
covariance update is 


Pk = (I'K k C k )P k (I-K k C k ) T + K k R k K k 

( 3 - 64 ) 

The selection of Kk is used to minimize the weighted sum of the diagonal elements of this 
error covariance matrix. Therefore, the cost function is 


J k =E[e k Se k ] 


( 3-65 ) 
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where S is any positive semidefinite matrix (i.e. I). Hence the cost function is just the trace 
of the error covariance matrix, which would be the same as minimizing the length of the 
estimation error vector. 


J k = trace[P k ] 


( 3-66 ) 


To determine the value of Kk that provides a minimum, it is necessary to create the 
Jacobian of the cost function with respect to the gain and set it equal to zero. Since it is 
known that [23] 


then 


a T 

— [trace(ABA )] = 2AB 
3 A 


(3-67) 


0 = -2(I-K k C k )P k C k +2K k R k 


( 3-68 ) 


Solving for Kk results in 


K k =P k ‘C k [C k P k X k +R k r' 


( 3-69 ) 


Gelb notes that the value of Kk calculated by using this equation is optimal and can be 
proven so through the examination of the Hessian of the cost function (i.e. Hessian of J k is 
positive semidefinite). 


( 3-70 ) 
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Using equation ( 3-69 ) and equation ( 3-64 ), the optimized value of the estimation error 
covariance matrix is calculated as 

p k = p k - p k C k[C k P k C][ + RJ‘‘C k P k 
= [I-K k CJP k 

(3-71) 

The state estimation and error covariance are extrapolated from one time step to another by 


*k “ Ak-l^k-1 


P k — ^k-l P k-l^k-l + Qk-I 


(3-72) 


It is helpful to see the discrete time Kalman filter variables in a graphical timing diagram. 
This helps to visually understand the steps needed in Kalman filtering. 


C k -1, Rk-1 Ck» Rk 



Figure 3-11. Discrete Kalman Filter Timing Diagram [23] 
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A summary of the discrete time Kalman Filter Equations can be found below. 



Kalman Filter Equations 

System Model 

x k = A k -. x k -. + w k -P w k _, ~N(0,Q k ) 

Measurement Model 

z k=C k x k +v k , v k ~ N(0,R k ) 

Initial Conditions 

E[x(0)]=x 0 , e[(x(0)-x 0 Xx(0)-x 0 )' ]=P 0 

Other Assumptions 

E[w k v J 'J = 0 forallj.k 

Time 

Update 

State Estimation Extrapolation 

X k = A k-l*k-l W k-I 

Error Covariance Extrapolation 

^k = ^k-l^k-l^k-l "*"Qk-l 

Measurement 

Update 

Kalman Gain Matrix 

K l =p.c;[c l p k c; + Rj' 

State Esitmate Update 

x k =x k +K k [z k -C k x k ] 

Error Covariance Update 

P,-[l-K,C,fc 


Table 3-4. Summary of Discrete Kalman Filter Equations [23] 


As Welsh and Bishop [26] describe it, the Kalman filter estimates a process by using a 
form of feedback control; the filter estimates the process state at some time and then 
obtains feedback in the form of noisy measurements. Therefore, the Kalman filter 
equations can be dividend into a two stage algorithm, a time update group and a 
measurement update group. The time update portion involves the forward projection of the 
current state and error covariance estimates to gain the a priori estimates needed for the 
measurement update. In the measurement update, the so called feedback occurs allowing 
for changes based on the new measurements. This new knowledge improves upon the a 
priori estimates and forms the improved a posteriori estimates. 
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Welsh and Bishop [26] liken this algorithm to a predictor (time update) corrector 
(measurement update) algorithm. Gelb has also named the two stages of the Kalman filter 
algorithm in a similar fashion as the extrapolation stage (time update) and the update stage 
(measurement update). The time update equation can be found in rows 5 and 6 of Table 
3-4 and the measurement update equations can be found in rows 7, 8, and 9 of Table 3-4. 


Notice that the time update equations project the state and error covariance estimates 
forward from time step k-1 to k, while the measurement update equations all work at time 
step k. In the measurement update the first step is to calculate a new Kalman gain Kk. The 
next step is to take a measurement of the process to get Zk. Then with this measurement, an 
a posteriori state estimate, x k , can be calculated. Then the final step of the Kalman filter 
algorithm iteration is to calculate the a posteriori error covariance matrix. The next 
iteration starts by using the last iteration’s a posteriori estimates as the new iteration’s a 
priori estimate. The recursive nature of the Kalman filter algorithm provides a large 
computational improvement on the Wiener filter, which is designed to operate on all of the 
data directly for each estimate [26]. 



Time Update 
(“Predict”) 


Measurement Update 
(“Correct”) 



Figure 3-12. Predictor - Corrector Model [26] 
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A block diagram of the plant and estimator can be found in Figure 3-13, which shows the 
State Estimate Extrapolation and the State Estimate Update as listed in Table 3-4. 


v 


k 


Measurement 



Figure 3-13. Block Diagram of Discrete System and State Estimator 


The Kalman gain, K k , in Figure 3-13 is calculated using the equations from Table 3-4 and 
is show in Figure 3-14. It is important to note that in the example used in this thesis, the 
deterministic input into the estimator is a zero vector since no disturbance measurements 
are available for state estimation. Figure 3-14 shows the Error Covariance Extrapolation, 
Error Covariance Update, and the Kalman Gain Computation as shown in Table 3-4. 
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Error 



Error Covariance Extrapolation 

Figure 3-14. Block Diagram Kalman Gain Computation 


The measurement noise covariance, Rk, is usually known since some measurements are 
taken. This allows for the calculation of Rk prior to the operation of the filter or at the 
beginning of the operation in some off-line process. On the other hand, the measurement 
of the process noise covariance is more difficult to calculate. This is reasonable due to the 
fact that there is no way to observe the process that is being estimating; if there was, there 
would be no need for estimation in the first place. However, acceptable results can result if 
one “injects” enough uncertainty into the process via the selection of Q [26]. This 
statement may seem vague, but it is important to understand that the process noise 
covariance is specific to each application. Therefore, performance will change for different 
values of Q. Common engineering techniques were used to find an initial value for Q (see 
Section 3.5), which remained constant throughout the testing process in order to 


standardize the results. 
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Albeit the limitations for Q, superior filter performance can still be achieved through tuning 
of Q and R. This usually involves an off-line process, involving a separate Kalman Filter, 
referred to as system identification. Increasing the process noise covariance effectively 
increases the bandwidth of the filter, which improves its tracking capabilities at the expense 
of more noise transmission [27]. If the value of Q is small it represents the belief that the 
Kalman filter model is a good representation of the plant. If the value of Q is large that 
represents our belief that the filter model is a poor representation and that trust in the 
measurement must be increased. In the special case where both Q and R are constants, 
both the estimation covariance and the Kalman gains are guaranteed to stabilize quickly 
and remain constant. 


3.5 Initial Kalman Filter Parameter Calculations 

Standard engineering methods were used in the determination of unknown initial 
measurement noise covariance, R, process noise covariance, Q, and error covariance, P. 
The initial Q matrix was set equal to a diagonal matrix representing a standard deviation of 
5% of the steady state amplitude, A ss , for the corresponding state. 


Q = diag(A ss *0.05) : 


(3-73) 


The steady state amplitudes can be calculated prior to filter operation through time domain 


simulations of the filter model. 
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The initial R matrix was set equal a diagonal matrix representing a standard deviation of 
2% of the steady state amplitude for the corresponding plant state. 

R = diag(A ss *0.02)' 

(3-74) 

The initial P matrix was set to equal a diagonal matrix, whose elements represent the 
square of three times the steady state amplitude for each state. 

P = diag(A ss *3) 2 

(3*75) 

Since there is no knowledge of the plant’s initial conditions, the initial conditions for the 
observer model were set to zero, even though the actual initial conditions applied to the 
plant were not. The initial conditions were set to equal the product of the steady state 
amplitude and a random number, rand, defined by a Gaussian distribution with a mean of 
zero and a standard deviation of 1. Using this distribution allowed for the initial conditions 
to be 180° out of phase. 


jc _ 0 _ env = A ss * rand 


(3-76) 
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4 Results 


Three different sets of tests were performed: 1) The first set of tests were conducted to 
verify the performance of disturbance modeling within the filter model, ISS Model 1, 
2) The second set of tests show the increased range of disturbance frequencies under which 
the Kalman filter is able to operate by using an expanded ISS disturbance model, ISS 
Model 2 and 3) The third set of tests included Monte Carlo analysis to determine 
robustness and investigate factors which have the most influence on the estimation error. 
These factors include differences between the plant and filter models (parameter variation), 
imbalance disturbance amplitude variation, and ISS disturbance frequency variation. 

4.1 Performance Measures 

Performance will be evaluated by using different metrics to include a measure of percent 
amplitude error in estimation, the duration of error, estimation error standard deviation, and 
time to convergence. 

4.1.1 Estimation Percent Amplitude Error 

Estimation percent amplitude error is defined as the ratio of the 2-norm of the error over 
the 2-norm of the actual state at steady state multiplied by 100% 


% amplitude error = 


pc 


. — i. *100% = 

ML 



2 - * 100 % 


( 4-1 ) 
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where x is the actual state, x is defined as the estimated state, and the error, x eiT , is defined 
as the difference between the actual state and the estimated state. 


4.1.2 Error Duration 

Corresponding to the percent amplitude error there is also an “error duration”, which 
measures the time that is spent within a certain percent amplitude error range with respect 
to the duration of the simulation. The percent amplitude estimation error, which is 
calculated at each time step, is collected into bins which divide the total range of percent 
amplitude estimation error (0 to 100%) equally (e.g., bin 1:0 to 5%, bin 2: 6 to 10%, bin 3: 
1 1 to 15%, etc). Dividing the number of occurrences in a given bin by the total number of 
simulation data point collected gives the percent of simulation time that resulted in an error 
within that bin’s range. See Figure 4-1 for an example of the error duration plot. 


Step Plot of Amplitude Error Percentage Duration 


This point is interpreted 
as -32% of all sim time 
has an error between 0% 
and 5% 


0 10 20 30 40 50 60 70 80 90 100 

Amplitude Error Percentage 

Figure 4-1. Example of Error Duration Performance Metric 
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While the percent amplitude error measures the error magnitude, the error duration 
provides insight to the severity of the estimation error. For instance a time history of brief 
but large error spikes may produce a large percent amplitude error, but looking at the error 
duration plot will show that the estimation error is relatively benign and may still be 
acceptable to slower controllers with appropriate robustness. 


4.1.3 Estimation Error Standard Deviation Envelope 

Yet another useful metric is the estimation error standard deviation envelope created by 
plotting the +/- square root of the error covariance, P, time history for a given state 
[27][28][29]. Plotting the amplitude error and the estimation envelope provides an 
indication of how often the error is outside of one standard deviation of the expected or 
predicted error values. From this plot, the percent of time the error spends outside of the 
estimation envelope can be calculated to assess the quality of the estimation. Figure 4-2 is 
an example plot of estimation error with the superimposed error standard deviation 
envelope for a case of excellent estimation, while an example of unacceptable estimation 
can be found in Figure 4-3. 


68 


State: X f Error Standard Deviation Envelope (Zoomed In) (% w/n ■ 100%) 



Figure 4-2. Example Plot of Estimation Error and Error Standard Deviation 
Envelope for Excellent Estimation (100% within bounds) 


State: x ( Error Standard Deviation Envelope (Zoomed In) (% w/n ■ 39.8512%) 



Figure 4-3. Example Plot of Estimation Error and Standard Deviation Envelope for 
Unacceptable Estimation (-40% within bounds) 
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4.1.4 Time to Convergence 

Time to convergence is defined as the time it takes for the estimated state to converge to 
the actual state. This can be determining by the time it takes for the error covariance to 
settle to a steady state value. For example in Figure 4-3, the square root of the error 
covariance settles sometime between ~3 seconds. As confirmed by the time history plot of 
the actual state and estimated state for the same case found in Figure 4-4, the state estimate 
converges to the actual state within 2.5 seconds. 


State: x ( (I s * 10 sec) (Actual vs Estimate) 



Figure 4-4. Time History Showing Convergence within 2.5 Seconds 
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4.2 ISS Disturbance Models Used for Testing 

The first ISS disturbance model, called “ISS Model 1” is the 4 state ISS disturbance model 
described by equations ( 3-6 ) though ( 3-9 ) in Section 3. 1.2.2. This model captures one 
translational ISS disturbance frequency (0)d S t) and one rotational disturbance frequency 

(tttdsr)- 


The second ISS disturbance model, called “ISS Model 2” is the 8 state ISS disturbance 
model described by equations ( 3-10 ) through ( 3-13 ) in Section 3. 1.2.3. ISS Model 2 
represents the ISS disturbance as sums of two sinusoidal disturbances, with two 
translational and two rotational disturbance frequencies of <jD dst i, 0J ds , 2 , and to^ri, to dsr 2 , 
respectively. 


4.3 ISS Model 1 Test: Performance 

The first set of tests attempt to answer two important questions: 1) Will disturbance 
modeling, of both the rotor disturbance and/or the ISS disturbance, inside the observer 
model, allow for estimation of absolute rotor states from relative measurement corrupted by 
sensor noise? 2) How much of an improvement is made over an observer with no 
disturbance modeling within the observer model? 

To answer these two questions, the same disturbances will be applied to: 1) an observer 
model with rotor and ISS disturbances modeling, 2) an observer model with only rotor 
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disturbance modeling, 3) an observer model with only ISS disturbance modeling, and 4) an 
observer model with no disturbance modeling. 

4.3.1 ISS Model 1 Test Set-Up 

Testing parameters will be chosen using two different ratios: 1) Frequency Ratio (FR) and 
2) Amplitude Ratio (AR). FR is the ratio of the ISS disturbance frequency, (0d S , to the rotor 
disturbance frequency, Wdr. as shown in equation ( 4-2 ). 


( 4 - 2 ) 

This ratio is used as a guide to identify worst case conditions. A worst case scenario for 
separating a single relative measurement into its components will occur when those 
components have similar frequency content. Therefore, after the rotor disturbance 
frequencies have been chosen, testing will occur such that the ISS disturbance frequency 
will be 90%, 100%, and 1 10% of the rotor disturbance frequency. 

The rotor disturbance frequencies were chosen to equal the peak mode for each of the 
following transfer functions from ISS disturbance inputs to relative measurements: 1) ISS 
disturbance in the x-axis to relative measurement in the x-axis (d xs to x re |), 2) ISS 
disturbance in the y-axis to relative measurement in the y-axis (d ys to y re i), 3) ISS 
disturbance about the x-axis to relative measurement about the x-axis (d<|> xs to <t) xre i)-, and 4) 
ISS disturbance about the y-axis to relative measurement about the y-axis (d^ to <j> yr ei)- 
See Figure 4-5 for Bode plots and peak frequencies used to determine testing frequencies. 
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The “dof of focus” will be determined by using the peak mode in a certain axis. For 
example, the dof of focus will be in the x-axis when a 0.399 Hz disturbance is used. 


Bode Plot: All ISS Disturbance Inputs 



Figure 4-5. Bode Plot Used to Determine Testing Frequencies 


AR is the output contribution ratio or the ratio of the relative measurement content due to 
the rotor disturbance, y re i_dr, and due to ISS disturbance, y re i_d S - 


AR = 


y rel _ ds 
y rel _ dr 


(4-3) 

Since the system is linear, superposition can be used. First, the relative motion is measured 
when only the rotor disturbances act on the plant. Then the relative motion is measured 
when only the ISS disturbances act on the plant. The amplitude of the ISS disturbances, F s , 
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are adjusted until the ratio of the plant output due to only ISS disturbance, is 10% and 
100% of the plant output due to only rotor disturbances. See Figure 4-6 below. Note that 
all other disturbance amplitudes are set to 1. See Table 4-1 for further explanation. 


Rotor 
Disturbance 


uisiuroance r? 

I *r 

| v v V \ 



°* 


ISS 

Disturbance 

Only 



(y rel y rel_dr + y rel_ds^ 

Figure 4-6. Amplitude Ratio Components 


y rel dr 


y relds 


Figure 4-7 shows a flow chart describing the method used for determining the ISS 
disturbance force amplitudes necessary to produce the desired ARs. 



Figure 4-7. Flow Chart to Determine F s Necessary for Desired Amplitude Ratios 
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With the rotor disturbance frequencies determined from Figure 4-5, the F s values were 
determined to get the desired AR. Table 4-1 contains the parameters used for testing 1SS 
Model 1. 


Dot of 
Focus 


Disturbance 
Frequencies | 

1 

ISS Force/Torque Amplitudes 

FR 

<*u (Hz) 

<* r (Hz) 

AR 

F,.(N) 

Fy. (N) 

T*. (Nm) 

T*. (Nm) 

X 

0.9 

0.359 

0.399 

0.1 

819.18 

l 

1 

1 

X 

0.9 

0.359 

0.399 

1 

8192.2 

l 

1 

1 

X 

1.0 

0.399 

0.399 

0.1 

9.9342 

l 

1 

1 

X 

1.0 

0.399 

0.399 

1 

99.316 

l 

1 

1 

X 

1.1 

0.439 

0.399 

0.1 

608.06 

l 

1 

1 

X 

1.1 

0.439 

0.399 

1 

6081.2 

l 

1 

1 


y 

0.9 

0.091 

0.101 

0.1 

1 

22.272 

1 

1 

y 

0.9 

0.091 

0.101 

1 

1 

222.73 

1 

1 

y 

1.0 

0.101 

0.101 

0.1 

1 

0.45137 

1 

1 

y 

1.0 

0.101 

0.101 

1 

1 

4.5132 

1 

1 

y 

1.1 

0.111 

0.101 

0.1 

1 

9.1608 

1 

1 

y 

1.1 

0.111 

0.101 

1 

1 

91.66 

1 

1 


<t>x 

0.9 

0.540 

0.600 

0.1 

1 

1 

118340 

1 

<t»x 

0.9 

0.540 

0.600 

1 

1 

1 

1183800 

1 


1.0 

0.600 

0.600 

0.1 

1 

1 

950.29 

1 

<t>x 

1.0 

0.600 

0.600 

1 

1 

1 

9508 

1 

<t»x 

1.1 

0.660 

0.600 

0.1 

1 

1 

152230 

1 

<!>x 

1.1 

0.660 

0.600 

1 

1 

1 

1522900 

1 


<t>y 

0.9 

0.711 

0.790 

0.1 

1 

1 

1 

30856 

0>y 

0.9 

0.711 

0.790 

1 

1 

1 

1 

308650 

4>y 

1.0 

0.790 

0.790 

0.1 

1 

1 

1 

181.8 

<t>y 

1.0 

0.790 

0.790 

1 

1 

1 

1 

1818.4 

Oy 

1.1 

0.869 

0.790 

0.1 

1 

1 

1 

39414 


1.1 

0.869 

0.790 

1 

1 

1 

1 

394250 


Table 4-1. Testing Parameters Determined for Desired FR and AR 
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4.3.2 ISS Model 1 Test Results 

As a summary of the data presented later in this section, Table 4-2 below helps to clearly 
show the conclusions that can be drawn with the examination of all test cases, which is: All 
disturbance modeling is necessary in order to perform good estimation. Table 4-2 shows 
the percent amplitude error in estimation for a case where all disturbances are modeled 
within the filter and a case where none of the disturbances are modeled, where each percent 
error value was given for the case were AR = FR = 1 for the given dof of focus. 


Dof of 
Focus 

Percent Amplitude Error 

All Disturbance 
Models 

No Disturbance 
Models 

X 

2.86% 

39.95% 

y 

0.18% 

30.33% 

<t>x 

0.43% 

91.47% 

<t>y 

1.45% 

59.67% 


Table 4-2. Percent Amplitude Error in Estimation between a Filter Model with All 
Disturbances Modeled and a Filter Model with No Disturbance Modeling 


Table 4-2 shows that without all disturbance modeling, estimation within acceptable error 
bounds is not possible. The rest of this section will go into further detail and give data for 
the different levels of disturbance modeling fidelities, but the conclusion that all 
disturbance modeling is necessary is still the same. 


A comparison of the results from using all disturbance models versus only rotor 


disturbance model for each test case can be found in Table 4-3. 
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Percent Am 

plitude Error 

Dof of 


Disturbance 

Frequencies 


ISS Force/Torque Amplitudes 

x » 

Yr 

♦xr 


Focus 

FR 

ob. (Hz) 

°*Jr (HZ) 

AR 

F,.(N) 

Fy.(N) 

T<m (Nm) 

T„. (Nm) 

All d 
Models 

No d„. 

All d 
Models 

No d,„ 

All d 
Models 

No d ( „ 

All d 
Models 

No d„. 

X 

0.9 

0.359 

0.399 

0.1 

819.18 

1 

1 

1 

5.42 

14.24 

3.94 

4.15 

1.02 

1.02 

1.45 

1.45 

X 

0.9 

0.359 

0.399 

1 

8192.2 

1 

1 

1 

11.32 

32.94 

3.26 

3.28 

1.24 

1.24 

1.80 

1.80 

X 

1.0 

0.399 

0.399 

0.1 

9.9342 

1 

1 

1 

3.47 

3.47 

4.01 

4.01 

0.98 

0.98 

1.38 

1.38 

X 

1.0 

0.399 

0.399 

1 

99.316 

1 

1 

1 

2.86 

2.98 

4.44 

4.44 

0.91 

0.91 

1.42 

1.42 

X 

1.1 

0.439 

0.399 

0.1 

608.06 

1 

1 

1 

4.14 

11.14 

4.02 

4.21 

1.01 

1.01 

1.37 

1.37 

X 

1.1 

0.439 

0.399 

1 

6081 .2 

1 

1 

1 

7.60 

44.43 

4.38 

4.47 

1.25 

1.25 

1.37 

1.37 

y 

0.9 

0.091 

0.101 

0.1 

1 

22.272 

1 

1 

1.10 

1.11 

0.47 

3.94 

1.10 

4.94 

1.26 

1.26 

y 

0.9 

0.091 

0.101 

1 

1 

222.73 

1 

1 

1.17 

1.17 

2.42 

23.52 

2.19 

5.20 

1.42 

1.42 

y 

1.0 

0.101 

0.101 

0.1 

1 

0.45137 

1 

1 

1.09 

1.09 

0.20 

0.24 

0.80 

0.81 

1.24 

1.24 

y 

1.0 

0.101 

0.101 

1 

1 

4.5132 

1 

1 

1.10 

1.10 

0.18 

0.63 

0.64 

0.87 

1.24 

1.24 

y 

1.1 

0.111 

0.101 

0.1 

1 

9.1608 

1 

1 

1.09 

1.09 

0.32 

2.06 

0.85 

2.00 

1.28 

1.28 

y 

1.1 

0.111 

0.101 

1 

1 

91.66 

1 

1 

1.13 

1.13 

1.51 

13.90 

1.22 

2.56 

1.59 

1.59 

Ox 

0.9 

0.540 

0.600 

0.1 

1 

1 

118340 

1 

2.46 

2.57 

0.50 

0.65 

0.66 

3.80 

1.40 

1.45 

Ox 

0.9 

0.540 

0.600 

1 

1 

1 

1183800 

1 

0.31 

0.32 

0.11 

0.13 

0.15 

0.78 

0.88 

0.94 

Ox 

1.0 

0.600 

0.600 

0.1 

1 

1 

950.29 

1 

6.27 

6.37 

3.90 

3.99 

1.09 

1.25 

1.25 

1.25 

Ox 

1.0 

0.600 

0.600 

1 

1 

1 

9508 

1 

1.80 

1.80 

0.34 

0.35 

0.43 

0.56 
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0.869 
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1 

1 

1 
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4.12 
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10.72 


Table 4-3. Comparison between Using All Disturbance Models and Only Rotor 


Disturbance Model: Test Results 


The percent amplitude error resulting from the use of all disturbance models is always less 
than or equal to the percent amplitude error resulting from using a disturbance model that 
does not include ISS disturbance. This is true for every rotor state and every test case. 
Test cases where the percent amplitude error is the same, or similar, for both disturbance 
modeling fidelities occurs only when the AR is small (i.e. rows where AR = 0.1) or if the 
frequency of the disturbances is not a peak mode frequency in that dof (i.e. all non 
highlighted results). This is reasonable since these two cases (small AR and non-peak 
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mode excitation) would cause a small response, therefore it would be expected that ISS 
disturbance modeling is not necessary for those cases. 


A comparison of the results from using all disturbance models versus only ISS disturbance 
model for each test case can be found in Table 4-4. 
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Ta 


ble 4-4. Comparison between Using All Disturbance Models and Only ISS 
Disturbance Model: Test Results 


The percent amplitude error resulting from the use of all disturbance models is always less 
than or equal to the percent amplitude error resulting from using a disturbance model that 
does not include rotor disturbance. This is true for every rotor state and every test case. By 
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comparing Table 4-3 to Table 4-4, it is evident that the rotor disturbance modeling is most 
important in reducing the percent amplitude error. 


A comparison of the results from using all disturbance models versus no disturbance 
models for each test case can be found in Table 4-5. 
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Table 4-5. Comparison between Using All Disturbance Models and Only ISS 


Disturbance Model: Test Results 


The percent amplitude error resulting from the use of all disturbance models is always less 
than or equal to the percent amplitude error resulting from using no disturbance modeling. 
This is true for every rotor state and every test case. This is expected because without 
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disturbance modeling, the filter model has no knowledge of the disturbances that act on the 
plant. 


The time history plots of one case (x r with AR = 1 and FR = 0.9) are used to highlight the 
need for full disturbance modeling within the filter model. The results for all different 
disturbance modeling fidelity levels can be found in Table 4-6. 


Percent Amplitude Error (x r ) 

All Disturbances 
Models 

Only Rotor 
Disturbance Model 

Only ISS 

Disturbance Model 

No Disturbance 
Models 

11 . 32 % 

32 . 94 % 

49 . 57 % 

56 . 74 % 


Table 4-6. Performance Verification Test Case for x r (AR = 1, FR = 0.9) 


It can be concluded from the results shown in Table 4-6 that both the rotor disturbance 
and ISS disturbance models are necessary for improving estimation capabilities. The rotor 
disturbance seems to have the highest effect on improving estimation capabilities. A 
possible explanation of this may be the facts that: 1) the rotor disturbance in the filter 
model is collocated with the relative measurement sensor, 2) due to the ‘nature’ or the plant 
being used, the relative measurement is mostly comprised of the rotor motion. In 
explanation of fact 2, the plant is stiff everywhere except between the shroud and the rotor 
(see Figure 2-2). Consequently, regardless of whether the disturbance is acting on the rotor 
or the outside mass of ISS flex model, the majority of the relative motion will come from 
the motion of the rotor. The frequency content of the sensor measurement, x re j, shows that 
the rotor motion, which is at a frequency of 0.399 Hz, is more important than the shroud 
motion, which is at a frequency of 0.356 Hz. 
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PSD of Relative Measurement: x |#| 



Frequency (Hz) 


Figure 4-8. Sensor Measurement x rc | Time History (left) and PSD (right) 


From this reasoning, one would expect that marginal estimation could be achieved when a 
rotor disturbance model exists, regardless of whether or not the ISS disturbance model is 
implemented, as long as sensor measurement content does not include a large component 
with the ISS disturbance frequency. The Kalman gains for the rotor disturbance model can 
adjust in order to compensate for the excess motion caused by the “unknown” disturbance 
source; the ISS disturbance is “unknown” to the filter since it is not modeled within the 
filter. With this said, the rotor-disturbance-only model can be improved upon with the 
modeling of the ISS disturbance, especially for the case where the relative measurement 
contains motion at the ISS disturbance frequency. Another compensation method would be 
to force the motion of the shroud to minimize the difference between the estimated relative 
measurement and the actual relative measurement. Error in shroud state estimation is of no 


concern since it is not an ABS controller input. 
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The following plots show for each level of disturbance modeling fidelity: 1) A comparison 
of the 1 st 10 seconds (left) and the last 10 seconds (right) of the actual and estimated value 
of x r (Figure 4-9 through Figure 4-12). 2) The actual amplitude error and the error bounds 
produced by the square root of the error covariance value (Figure 4-13 though Figure 4-16). 


State: x f (1 st 20 sec) (Actual vs Estimate) 




Figure 4-9. Results of Implementation of All Disturbance Models 


State: x, (I* 20 tec) (Actual va Estimate) , State: x (Last 20 tec) (Actual va Estimate) 

1 x 10 r 




Figure 4-10. Results of Implementation of Only Rotor Disturbance Model 
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State: x (1 st 20 sec) (Actual vs Estimate) 3 State: x (Last 20 sec) (Actual vs Estimate) 

■ x 10 r 




Figure 4-11. Results of Implementation of Only ISS Disturbance Model 


State: x (I s1 20 sec) (Actual vs Estimate) 



State: x (Last 20 sec) (Actual vs Estimate) 

x 10 1 



Figure 4-12. Results of Implementation of No Disturbance Models 


The Kalman filter estimate converges very quickly (within 5 seconds) to the actual state in 
all four different disturbance modeling fidelity levels. As expected, when all of the 
disturbances that are applied to the plant are being modeled within the filter model, good 
estimation is achieved (Figure 4-9). Also, the estimation error is well within the one 
standard deviation envelope 100% of the time (Figure 4-13). 
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The lack of ISS disturbance modeling in the Only Rotor Disturbance Model (Figure 4-10) 
creates amplitude error between the actual and the estimated rotor state. The estimation 
error cycles in and out of being within one standard deviation of expected error (Figure 
4-14). This shows that there are some problems with the filter process when no ISS 
disturbance is being modeled within the filter model. 

When the rotor disturbance is not modeled within the filter model, there seems to be some 
phase error due to the lack of information of the rotor disturbance within the filter model 
(Figure 4-11). This causes the cycling estimation error (Figure 4-15) to stray further from 
the error standard deviation envelope and for longer periods of time when compared to the 
result for the Only Rotor Disturbance Model results shown in Figure 4-14. 

The No Disturbance Models test (Figure 4-12) resulted in both phase error and amplitude 
error, which combined to create the largest error. Poor estimation as evidenced by the 
increased occurrence of the estimation error exceeding the standard deviation envelope 
(Figure 4-16). 
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State: x ( Error Standard Deviation Envelope (Zoomed In) (% win = 100%) 



Figure 4-13. Standard Deviation Envelope Resulting from Implementation of 
All Disturbance Models 


State: x ( Error Standard Deviation Envelope (Zoomed In) (% w/n ■ 66.9407%) 



Figure 4-14. Standard Deviation Envelope Resulting from Implementation of 
Only Rotor Disturbance Model 


State: x ( Error Standard Deviation Envelope (Zoomed In) (% w/n = 39.8512%) 



Figure 4-15. Standard Deviation Envelope from Implementation of Only ISS 
Disturbance Model 


State: x ( Error Standard Deviation Envelope (Zoomed In) (% w/n * 37.0593%) 



Figure 4-16 Standard Deviation Envelope from Implementation of No 
Disturbance Models 
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The error duration plots continue to prove that without the use of all disturbance models 
within the filter, poor estimation will result. 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-17. Duration of Error of x r Estimation from Implementation of 
All Disturbance Models 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-18. Duration of Error of x r Estimation from Implementation of 
Only Rotor Disturbance Model 


Percentage of Sim Time ™ Percentage of Sim Time 


Step Plot of Amplitude Error Percentage Duration 



-19. Duration of Error of x r Estimation from Implementation of 
Only ISS Disturbance Model 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-20. Duration of Error of x r Estimation from Implementation of 
No Disturbance Models 
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As the disturbance modeling fidelity decreases, the percent amplitude error durations 
become more distributed. For instance, when all disturbances are modeled, -55% of the 
simulation time had an error of under 10%, while only -14% of the simulation time was 
under 10% error when using no disturbance models. 


4.3.3 ISS Model 2 Testing: Increased Disturbance Frequency 
Range Test 

These tests attempt to answer the question: Will the use of additional disturbance states 
increase the range of ISS disturbance frequencies under which the Kalman filter is able to 
operate? This question is answered by comparing the performance of ISS Model 1 to ISS 
Model 2. The rotor disturbance frequency was set to 0.399 Hz, when focusing on x and <|>y 
dofs, and set to 0.101 Hz when focusing on y and <t> x dofs. These frequencies are the peak 
modes in the following transfer functions: ISS disturbance in the x-axis to x re i measurement 
(d xs to x re i) and ISS disturbance in the y-axis to y re | measurement (d ys to y re i), respectively. 
See Figure 4-21 for more details. Since the peak frequencies are being used, the coupled 
motions (x and <|> y ) and (y and <|>x) will have the greatest amplitude when excited by those 
frequencies with respect to the other degrees of freedom. Each axis will be examined 
separately to help provide clear result from which sound conclusions can be drawn. 
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4.4 ISS Model 2 Test Set-Up 

The ISS translation and rotational disturbance frequencies for each test can be found in 
Table 4-7. In summary, the ISS translational disturbance frequencies were set to 90% of 
the rotor disturbance frequency for Test A and 1 10% of the rotor disturbance frequency for 
test B. 


Bode Plot: Translational TF with ISS Inputs 



Figure 4-21. Bode Plots Used For Rotor Disturbance Frequencies 


As for the rotational ISS disturbance frequencies, when focusing on x and <j> y , the rotational 
ISS disturbance frequency was centered on the peak mode frequency (0.79 Hz) of the 
transfer function from the rotational ISS disturbance about the y axis to the <t> yre i 
measurement (d<t, ys to <f) yr et)- The rotational ISS disturbance frequency was set to 90% of 
this peak value (0.71 1 Hz) for Test A and 1 10% of this peak value (0.869 Hz) for Test B. 
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When focusing on y and (|) x , the rotational ISS disturbance frequency was centered on the 
peak mode frequency (0.6 Hz) of the transfer function from the rotational ISS disturbance 
about the y axis to the 0 yre i measurement (d 0 ys to <J) yre i)- The rotational ISS disturbance 
frequency was set at 90% of this peak value (0.54 Hz) for Test A and 1 10% of this peak 
value (0.66 Hz) for Test B. 


Bode Plot: Rotational TFs with ISS Disturbance Inputs 



Figure 4-22. Bode Plots Used For ISS Rotational Disturbance Center Frequencies 


ISS Model 2 has both “90%” and “110%” frequencies modeled, while ISS Model 1 has 
only the “90%” frequencies modeled. All of these frequencies are summarized in Table 
4-7. A second set of tests were run incorporating a rotor spin-up from 0 to 0.7 Hz over 300 


seconds. 


91 


Units: Hz 

Rotor Spin 
Frequency 

(C0.pln) 

ISS disturbance 
Frequencies Applied 
to Plant 

ISS Model 1 

(4 d,ss Filter 
States) 

ISS Model 2 

(8 diss Filter States) 

COt 

applied 

(Or 

applied 

90% Frequencies 

Both 90% and 110% Frequencies 

^dstl 

^dsrl 

w dst1 

W dst2 

® dsrl 

®dsr2 

x and <)>y 

Test A (90% Freq) 

0.399 



0.359 

0.711 

0.359 

0.439 

0.711 

0.869 

Test B (110% Freq) 

0.439 

0.869 

y and <)>„ 

Test A (90% Freq) 

0.101 

0.091 

0.54 

0.091 

0.54 

0.091 

0.111 

0.54 

0.66 

Test B (110% Freq) 

0.111 

0.66 

x and 4> y 

(w/spln-up) 

Test A (90% Freq) 

0.399 

0.359 

0.711 

0.359 

0.711 

0.359 

0.439 

0.711 

0.869 

Test B (110% Freq) 

0.439 

0.869 

y and 4> x 

(w/spin-up) 

Test A (90% Freq) 

0.101 

0.091 

0.54 

0.091 

0.54 

0.091 

0.111 

0.54 

0.66 

Test B (110% Freq) 

0.111 

0.66 


Table 4-7. Test Matrix 


To interpret of Table 4-7, the testing focusing on x and <J> y proceeded as follows: 

1) Focusing on x and <J) y , Test A was conducted by exciting the 32 state plant with 
4 rotor disturbances with a frequency equal to the spin rate (0.399 Hz), and 4 
ISS disturbances. The translational 1SS disturbances are pulse trains with a 
frequency of 0.359 Hz while the rotational ISS disturbances are pulse trains 
with a frequency of 0.711 Hz. These two frequencies have been labeled the 
“90%” frequencies. The relative measurement, x re i, is fed into two different 
Kalman Filters. The first Kalman filter includes the rotor disturbance model 
and the ISS (disturbance) Model 1 within the filter dynamics. ISS Model 1 
contains information about the two ISS disturbance frequencies applied (0.359 
and 0.711 Hz, the 90% frequencies). The second Kalman Filter includes the 
rotor disturbance model and the ISS (disturbance) Model 2 within the filter 
dynamics. ISS Model 2 contains knowledge of the same two ISS disturbance 
frequencies as modeled in ISS Model 1, the 90% frequencies, but also contains 
information about 2 additional ISS disturbance frequencies (0.439 and 0.869 
Hz, the 110% frequencies). State estimation error, x^ for ISS Model 1 and 
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x 8err f° r ISS Model 2, are computed by taking the difference from the estimated 
states from ISS Model 1, x 4 , and the estimated states from ISS Model 2, x g , 
and actual state, x, in order to determine performance. 

2) Again focusing on x and <|>y. Test B was conducted by exciting the 32 state plant 
with 4 rotor disturbances with a frequency equal to the spin rate (0.399 Hz), and 
4 ISS disturbances as was done in Test A, but this time, the translational ISS 
disturbance pulse trains are input at a frequency of 0.439 Hz while the rotational 
ISS disturbances are input at a frequency of 0.869 Hz. These two frequencies 
have been named the “1 10%” frequencies. The output of the plant, the relative 
measurement (x re i), is fed into the same two Kalman Filters as in Test A, and the 
process used to determine performance is also the same 

These same two tests are also run while focusing on y and <f> x dofs. The testing algorithm is 
show in Figure 4-23. 
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Figure 4-23. Testing Algorithm for Improving Operational Bandwidth 


The expectation for Test A is that both Kalman filters should produce similar errors. ISS 
Model 2 may have slightly higher error due to the fact that it is modeling a disturbance 
frequency that does not exist, but the filter gains that error state so the amplitude of the 
force with the unapplied frequency is very small. For Test B, the expectation is that the 
error in ISS Model 2 should be much less than in ISS Model 1. This is due to the fact that 
ISS Model 1 does not have the “110%” ISS disturbance frequency information modeled 
within the filter dynamics, while the ISS Model 2 does. See Figure 4-24 for a logic flow 
diagram that sums up the previous discussion. 
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Figure 4-24. Logic Flow Diagram for Testing Regimen 


The test results found in Section 4.4.1 confirm this hypothesis. 


4.4.1 ISS Model 2 Test Results: Focusing on x and <t> y 

When focusing on x and <j) y , the following frequencies were used for Test A and Test B. 


Test A(o d „ 


Test B a) d „ 
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(0*2=0.439 Hz 

. Mot: ISS Disturbi 

( 0 ,, =0.711 Hz 
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Fraquancy (H2) 

--.J 

( 0 ,,= 0.711 Hz 


Figure 4-25. Disturbance Frequencies Used when Focusing on x and (j^ 


The results from focusing on coupled x and <J> y rotor motion, without spin-up, can be found 


in Table 4-8. 
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Axis of Focus : 

x and ( 

h 

Test: 

A 

B 

ISS Model #: 

1 

2 

1 

2 

% imp 

Rotor 

States 

% Amp Error (x r ) 

6.57% 

6.85% 

17.18% 

3.62% 

78.93% 

% Amp Error (y r ) 

1.52% 

1.53% 

1.20% 

1.18% 

2.16% 

% Amp Error (<j>x r ) 

0.74% 

0.74% 

1.08% 

1.03% 

5.01% 

% Amp Error (<j>y r ) 

1.21% 

1.13% 

7.02% 

0.67% 

90.44% 


Table 4-8. Test Results Focusing on x and <)>y (no spin-up) 


As predicted, the percent amplitude error, calculated with equation (4-1 ), is nearly 
identical for both ISS Model 1 and ISS Model 2 for the Test A case. Also as expected, 
during Test B, ISS Model 2 shows a percent amplitude error improvement over ISS Model 
1 for all rotor states, and the largest improvements are for the x r and <t> yr dofs as predicted. 
The percent improvement (% imp) is calculated as the difference between the ISS Model 1 
result and the ISS Model 2 result divided by the ISS Model 1 result, therefore it is a 
measure of percent improvement in percent amplitude error over the ISS Model 1 error. 
This shows that the range of ISS disturbance frequencies that the filter will operate under 
can be increased by expanding and improving the disturbance model within the filter 
dynamics. 

As an example of the estimation improvement made by using ISS Model 2, the estimation 
of x r during Test B will be compared between ISS Models 1 and 2. The time history plot 
comparisons are shown in Figure 4-26 and Figure 4-27, the estimation error envelope 
comparisons are shown in Figure 4-28 and Figure 4-29, and the error duration comparisons 
are shown in Figure 4-30 and Figure 4-31 
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Figure 4-26. Time History of Actual, Estimated, and Error for x r Using ISS Model 1 
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Figure 4-27. Time History of Actual, Estimated, and Error for x r Using ISS Model 2 
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Both observers converge very quickly, but the steady state error for ISS Model 2, shown in 
Figure 4-27, is much smaller than when using ISS Model I, shown in Figure 4-26. 


State: x Standard Deviation Envelope (Zoomed In) (% w/n ■ 40.4472%) 



Figure 4-28. Estimation Error and Error Standard Deviation Envelope of x r 
Estimation with ISS Model 1 


State: X ( Standard Deviation Envelope (Zoomed In) {% w/n - 99.99%) 



Figure 4-29. Estimation Error and Error Standard Deviation Envelope of x r 
Estimation with ISS Model 2 


The estimation error for ISS Model 1, shown in Figure 4-28, stays within the standard 
deviation envelope only 40.44% of the time while the estimation error for ISS Model 2, 
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shown in Figure 4-29, stays within the standard deviation envelope nearly 100% of the 
time, showing the superior estimation capabilities when using ISS Model 2. 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-30. Duration of Error of x r Estimation with ISS Model 1 


Step Plot of Amplitude Error Percentage Duration 



Amplitude Error Percentage 


Figure 4-31. Duration of Error of x r Estimation with ISS Model 2 
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Figure 4-31 shows that state estimation is close at all times when using ISS Model 2. In 
fact over 75% of the simulation time has an estimation amplitude error less than 5%, and 
nearly 100% of the simulation time experiences an error of less than 20%. However, when 
using ISS Model 1, Figure 4-30, there is a greater distribution of error durations. In fact, a 
significant amount of the simulation time has a resulting estimation amplitude error of 
greater than 30%. 

Results with a rotor spin-up can be found in Table 4-9 below. 



Axis of Focus : 

x and (|) y w/ spin-up 

Test: 

A 

B 

ISS Model #: 

1 

2 

1 

2 

% imp 

Rotor 

States 

% Amp Error (x r ) 

20.20% 

20.89% 

50.53% 

9.14% 

81.90% 

% Amp Error (y r ) 

0.01% 

0.01% 

0.05% 

0.01% 

89.42% 

% Amp Error (fa r ) 

0.00% 

0.00% 

0.00% 

0.00% 

20.00% 

% Amp Error (<^ r ) || 2.84% 

2.75% 

7.30% 

1.57% 

78.45% 


Table 4-9. Results Focusing on x and (with spin-up) 


The hypothesis also holds true during rotor spin-up from 0 to 0.7 Hz over 300 seconds. In 
fact, for rotor state x r during Test B a substantial improvement is made from an amplitude 
error of 50.53% to an amplitude error of 9.14% with the use of ISS Model 2. The large 
percent improvement values for y r and <)) xr , should be disregarded as it is a numerical 
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artifact, meaning that the values are very small and therefore any small changes will 
produce a large improvement, even though the actual improvement is minuscule. 

As an example of the estimation improvement made by using ISS Model 2 with time- 
varying inputs and dynamics, the estimation of x r during Test B will be compared between 
ISS Models 1 and 2. The time history plot comparisons are shown in Figure 4-32 and 
Figure 4-33, the estimation error envelope comparisons are shown in Figure 4-34 and 
Figure 4-35, and the error duration comparisons are shown in Figure 4-36 and Figure 4-37. 


Displacement (m) Displacement (m) v Displacement (m) Displacement (m) 
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State: x ( (1** 20 sec) (Actual va Estimate) 



, State: x (Last 20 sec) (Actual vs Estimate) 
x 10 i 



State: x ( (I* 1 20 sec) (Actual vs Error) 
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Figure 4-33. Time History of Actual, Estimated, and Error for x r Using ISS Model 2 
(with rotor spin-up) 
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Both observers converge very quickly, but the steady state error for ISS Model 2, Figure 
4-33, is much smaller than when using ISS Model 1, Figure 4-32. The error when using 
ISS Model 1 is mainly due to amplitude estimation error rather than phase estimation error. 


State: x r Standard Deviation Envelope (Zoomed In) (% w In - 42.6011%) 



Figure 4-34. Estimation Error and Error Standard Deviation Envelope of x r 
Estimation with ISS Model 1 (with rotor spin-up) 


State: x f Standard Deviation Envelope (Zoomed In) (% wfn ■ 99.996%) 



Figure 4-35. Estimation Error and Error Standard Deviation Envelope of x r 
Estimation with ISS Model 2 (with rotor spin-up) 
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The estimation error for ISS Model 1, Figure 4-34, stays within the standard deviation 
envelope only 42.5% of the time while the estimation error for ISS Model 2, Figure 4-35, 
stays within the standard deviation envelope nearly 100% of the time, showing the superior 
estimation capabilities when using ISS Model 2. 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-36. Duration of Error of x r Estimation with ISS Model 1 (with 
rotor spin-up) 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-37. Duration of Error of x r Estimation with ISS Model 2 (with 
rotor spin-up) 
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Figure 4-37 shows that state estimation is close at all times when using ISS Model 2. In 
fact nearly 70% of the simulation time has an estimation amplitude error less than 5%, and 
nearly 100% of the simulation time experiences an error of less than 20%. However, when 
using ISS Model 1, Figure 4-38, there is a greater distribution of error durations. In fact, a 
significant amount of the simulation time has a resulting estimation amplitude error of 
greater than 30% 


4.4.2 ISS Model 2 Test Results: Focusing on y and <J> X 

When focusing on y and <(> x , the following frequencies were used for Test A and Test B. 



Figure 4-38. Disturbance Frequencies used when focusing on y and <t>x 
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The results from focusing on coupled y and rotor motion, without spin-up, can be found 
in Table 4-10. 



Axis of Focus : 

y and <|> x 

Test: 

A 

B 

ISS Model #: 

1 

2 

1 

2 

% imp 

Rotor 

States 

% Amp Error (x r ) 

2.79% 

2.79% 

2.01% 

2.00% 

0.24% 

% Amp Error (y r ) 

0.61% 

0.75% 

4.72% 

0.78% 

83.38% 

% Amp Error (<t>xr) 

0.58% 

0.58% 

3.67% 

0.44% 

87.99% 

% Amp Error (<t»y r ) 

0.91% 

0.92% 

0.60% 

0.59% 

0.41% 


Table 4-10. Results Focusing on y and (no spin-up) 


Again as predicted, the percent amplitude error is nearly identical for the Test A case for 
both ISS Model 1 and ISS Model 2. Also, the Test B shows an improvement, in the 
percent amplitude error for all rotor states, and shows large improvements for y r and <t> xr as 
predicted. 

As an example of the estimation improvement made by using ISS Model 2, the estimation 
of y r during Test B will be compared between ISS Models 1 and 2. The time history plot 
comparisons are shown in Figure 4-39 and Figure 4-40, the estimation error envelop 
comparisons are shown in Figure 4-41 and Figure 4-42, and the error duration comparisons 
are shown in Figure 4-43 and Figure 4-44. 


Displacement (m) Displacement (m) 
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State: y f (I - SO sec) (Actual vs Estimate) 



State : y f (Last SO sec) (Actual vs Estimate) 



State: y f (1 st 60 sec) (Actual vs Error) 



State: y f (Last SO sec) (Actual vs Error) 



Figure 4-39. Time History of Actual, Estimated, and Error for y r Using ISS Model 1 
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Figure 4-40. Time History of Actual, Estimated, and Error for y r Using ISS Model 2 
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Event though the estimation using ISS Model 1 is good. Figure 4-40 shows that the steady 
state error can be nearly eliminated with the use of increased disturbance modeling found 
in ISS Model 2. 


State: y ( Standard Deviation Envelope (Zoomed In) (% w In - 87.2483%) 
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Figure 4-41. Estimation Error and Error Standard Deviation Envelope of 
y r Estimation with ISS Model 1 


State: y f Standard Deviation Envelope (Zoomed In) (% w/n ■ 100%) 



Figure 4-42. Estimation Error and Error Standard Deviation Envelope of 
y r Estimation with ISS Model 2 
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Improved estimation using ISS Model 1 is also evident from the fact that estimation error 
does not stray outside the bounds of the estimation error standard deviation envelope 
(Figure 4-42). 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-43. Duration of Error of y r Estimation with ISS Model 1 


Step Plot of Amplitude Error Percentage Duration 



Figure 4-44. Duration of Error of y r Estimation with ISS Model 2 


109 


Figure 4-44 shows that the occurrence of high estimation error can be decreased by using 
the expanded disturbance modeling found in ISS Model 2. 



Axis of Focus : 

y and <|> x w/ spin-up 

Test: 

A 

B 

ISS Model #: 

1 

2 

1 

2 

% imp 

Rotor 

States 

% Amp Error (x r ) 

0.00% 

0.00% 

0.04% 

0.04% 

0.00% 

% Amp Error (y r ) 

0.81% 

0.93% 

4.58% 

0.81% 

82.34% 

% Amp Error (fa r ) 

0.82% 

0.80% 

2.66% 

0.71% 

73.32% 

% Amp Error (<t>y r ) 

0.00% 

0.00% 

0.00% 

0.00% 

0.00% 


Table 4-11. Results Focusing on y and <t>x (with spin-up) 


The hypothesis also holds true during the rotor spin-up from 0 to 0.7 Hz over 300 seconds. 
The results are shown in Table 4-11. 


4.5 Monte Carlo Analysis: Robustness Test Set-Up 

In order to test robustness to parameter uncertainty, each plant stiffness value was 
independently allowed to deviate from its nominal value where the deviations were defined 
by Gaussian distributions with a mean of zero and a 3a value of 20% of the nominal value. 
To test robustness to rotor imbalance disturbance amplitude, the components which make 
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up the rotor imbalance disturbance (M, e, a), see Figure 3-2, were allowed to individually 
deviate from their nominal value where the deviations were defined by Gaussian 
distributions with a mean of zero and a 3a value of 2 kg, 0.1 m, and 0.05 radians for M, e, 
and a, respectively. To test robustness to ISS disturbance frequency uncertainty, the 
disturbance frequencies were allowed to deviate from their nominal values where the 
deviations were defined by Gaussian distributions with a mean of zero and a 3a value of 
20% of the nominal value. 

Special care was taken in defining the distributions of the different deviations in order to 
prevent impossible deviations such as negative rodent mass, M. The Monte Carlo testing 
involved 1000 test runs of each uncertainty category using ISS Model 2. Other than the 
allowed deviations on the plant, nothing else in the simulation or the observer model was 
changed. The nominal spin frequency used equals 0.7 Hz, while the nominal ISS 
translational and rotational disturbance frequencies were chosen as the peak mode 
frequencies from the following two transfer functions, respectively: 1) d ys to y re i and 2) d^s 
to <t>xrei (see Figure 4-38). 


4.5.1 Monte Carlo Analysis: Parameter Uncertainty Results 

The following distributions are the result of the variation on the plant stiffness parameters 


as defined in Section 4.5. 


Ill 


K w< Nominal -928 Nominal » 928 Nominal -7630 Nominal - 7830 




Figure 4-45. Distribution of Translational and Rotational Stiffness between Rotor and 
Shroud (top row) and Shroud and ISS Mass 1 (bottom row) 



Figure 4-46. Distribution of Translational and Rotational Stiffness between ISS Mass 
1 and ISS Mass 2 (top row) and ISS Mass 2 and Inertial (bottom row) 
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The variation of stiffness values forces a distribution of the plant’s frequency response. 
The modal frequency distributions for the direct transfer functions (i.e. d xs to x re i) caused by 
the change in the stiffness parameters of the plant are shown in Figure 4-47 through Figure 
4-50. The nominal Bode plots have been superimposed on the modal frequency 
distributions to show the variation from the nominal plant modal frequencies caused by the 
distribution of the plant stiffness parameters. Figure 4-47 through Figure 4-50 also 
provide information on which frequencies are more likely to be effected; that is, which 
frequencies will shift or appear/disappear due to variations in the stiffness parameters. 


TF Bode and Distribution of Modal o» Due to aK 

x n 



Figure 4-47. Modal Frequency Distribution Caused by AK (for d xs to x rc |) 
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Figure 4-48 and Figure 4-49 show that the higher modal frequencies will have a larger 
standard deviation from the nominal value than the lower frequencies. 

TF Bode and Distribution of Modal <» DuetoAK 



Figure 4-48. Modal Frequency Distribution Caused by AK (for d ys to y rc i) 


TF, Bode and Distribution of Modal &> Due to AK 

(pX n 



Figure 4-49. Modal Frequency Distribution Caused by AK (for d,(, xs to <J>xrei) 
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TF, Bode and Distribution of Modal «> Due to aK 

♦y n 



Figure 4-50. Modal Frequency Distribution Caused by AK (for d^ to <j>y rc i) 


Figure 4-49 shows that the high frequencies are likely to shift away from the nominal value 
since the mean is not equal to the nominal value. Also in Figure 4-50, a new modal 
frequency appears around 0.58 Hz due to variations in plant stiffness parameters. 

The Monte Carlo analysis results, testing estimation performance under conditions of 
parameter uncertainty can be found in Figure 4-51 below. 
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% Amplitude Error (y |o|{M ) ^ - 6.3 / ti • 3.8 



% Amplitude Error 




Figure 4-51. Monte Carlo Results Testing Parameter Uncertainty (no spin-up) 


The largest percent amplitude error occurs in the translation of the rotor in the x axis. The 
estimation of the other dofs are very accurate. 

It is important to note that the process noise covariance matrix, Q, does not change from 
the nominal Q during the Monte Carlo simulations. Since Q is determined by the level of 
uncertainty in the observer model, by increasing Q as the difference between the plant 
model and the observer model increased due to parameter changes, improved estimation 
would be expected. Time history plots, displayed in Figure 4-52 through Figure 4-55, were 
created using a set of parameters that produced the approximate mean amplitude error 
values, resulting from the Monte Carlo analysis, for all four motions of the rotor (x r , y r 
translations and 0 xr , 0 yr rotations). These parameter values produced a percent amplitude 
error of x r — * 1 1.79%, y r — * 6.3%, 0 xr —* 2.5%, and 0 yr — ► 1.7%. 


Displacement (m) Displacement (m) 
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Figure 4-52. x r Estimation 
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Figure 4-53. y r 
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Figure 4-54. ((^r Estimation 
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Figure 4-55. <j>y r Estimation 


It is clear from Figure 4-52 through Figure 4-55 that the error is caused purely by 
amplitude differences between the actual and estimated states and not by phase lag or lead. 


Occurrences out of 1000 Occurrences out of 1000 


118 


The Monte Carlo analysis results, with the time-varying case of rotor spin-up, can be found 
in Figure 4-56 below. 


% Amplitude Error (x (oto| ) 11.3 / o- 9.1 % Amplitude Error (y |o|w ) ^ - 6.4 / a - 4 





% Amplitude Error 



Figure 4-56. Monte Carlo Results Testing Parameter Uncertainty (with spin-up) 


By comparing Figure 4-51 with Figure 4-56, it is evident that the Kalman filter 
performance is nearly identical regardless of time varying dynamics or time varying inputs. 
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4.5.2 Monte Carlo Analysis: Rotor Disturbance Amplitude 
Uncertainty Results 

The following distributions are the result of the variation of the parameters that determine 
the amplitude of the rotor disturbance as defined in Section 4.5. 



Figure 4-57. Distribution of Rotor Disturbance Amplitude Parameters 


The Monte Carlo analysis results, testing estimation performance under conditions of rotor 
disturbance amplitude uncertainty can be found in Figure 4-58. 



% Amplitude Error 



0 i 

0.214 0.2145 0.215 0.2155 0.216 0.2165 0.217 

% Amplitude Error 



Figure 4-58. Monte Carlo Results Testing Rotor Disturbance Amplitude 
Uncertainty (no spin-up) 
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5 Conclusions 

This thesis develops a state estimation algorithm for the Centrifuge Rotor (CR) system 
where only relative measurements are available with limited knowledge of both rotor 
imbalance disturbances and ISS thruster disturbances. A Kalman filter is applied to a plant 
model augmented with sinusoidal disturbance states used to model the effect of the ISS 
thrusters on the CR relative motion measurement. The sinusoidal disturbance states 
compensate for the lack of the availability of plant inputs for use in the Kalman filter. 
Testing confirms that complete disturbance modeling is necessary to ensure reliable 
estimation. Further testing goes on to show that increased estimator operational bandwidth 
can be achieved through the expansion of the disturbance model within the filter dynamics. 
In addition, Monte Carlo analysis shows the varying levels of robustness against defined 
plant/filter uncertainty variations. 

Chapter 2 provided a problem overview and a concise description of the CR system. This 
included a description of the simplified model as well as the list of assumptions necessary 
for simplification. The model used for analysis and testing included a 4 mass (Rotor, 
shroud, ISS Mass 1, and ISS Mass 2), 16 degree of freedom, time-varying system. From 
this model, linearized equations of motion were derived. 

In Chapter 3 a detailed description of how the disturbance models were implemented 
within the filter dynamics was formulated. This formulation required the derivation of both 
rotor and ISS disturbances. The rotor disturbance was derived as a function of imbalance 


geometry, mass/inertia, and spin rate, while the ISS disturbance modeling was done though 
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% Amplitude Error 


% Amplitude Error (♦ MIOtol ) n * 3.2 1 o ■ 2.6 



% Amplitude Error 


% Amplitude Error (y rotor ) M-* 14.7 / o - 11.7 




% Amplitude Error 


Figure 4-64. Monte Carlo Results Testing Combination of All Uncertainties 
(with spin up) 


Relatively large errors result from the combination of all three uncertainty categories. 
Error in estimating <]) xr and <J> yr rotations remain relatively low while the errors in both x r and 
y r translation are relatively high. The two uncertainties having the most effect on 
estimation error are the plant parameter uncertainty and the ISS disturbance frequency 
uncertainty, both of which create frequency disparities between the plant and the filter 


models. 
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4.5.4 Monte Carlo Analysis: Combination of All Uncertainties 
Results 

Monte Carlo analysis was conducted combining the plant parameter uncertainty, the rotor 
disturbance amplitude uncertainty, and the ISS disturbance frequency uncertainty in order 
to determine the error for a case with all uncertainties acting at the same time. The results 
for this analysis can be found in Figure 4-63 for the time-invariant case, and in Figure 4-64 
for the rotor spin-up case. 


% Amplitude Error (x ) ji - 23.6 / a ■ 22.3 
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% Amplitude Error (y, 0t0f ) |i" 14.3 /a" 11.5 



% Amplitude Error 



% Amplitude Error 
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Figure 4-63. Monte Carlo Results Testing Combination of All Uncertainties (no 
spin up) 
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Figure 4-61 shows that the largest errors are caused by ISS disturbance frequency 
uncertainly. Therefore, it will be important to determine the ISS disturbance frequencies of 
concern and to expand the ISS disturbance model to capture all of them. 


The Monte Carlo analysis result, considering a time-varying plant and disturbance inputs, 
can be found in Figure 4-62. 


% Amplitude Error ( x rot0f ) n" 7.13 /o" 3.74 
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% Amplitude Error (♦ yTOto( ) n" 1.41 / a- 0.02 



Figure 4-62. Monte Carlo Results Testing ISS Disturbance Frequency 
Uncertainty (with spin up) 


Again, similar results are found between the time-invariant and time-varying Monte Carlo 
analysis, supporting the conclusion that the Kalman Filter operates similarly for both cases. 
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4.5.3 Monte Carlo Analysis: ISS Disturbance Frequency 
Uncertainty Results 

The distributions, shown in Figure 4-60, are the result of the variation of both the 


translational and rotational ISS disturbance frequencies as defined in Section 4.5. 



Figure 4*60. Distribution of Translational and Rotational ISS Disturbance 
Frequencies 


The Monte Carlo analysis results, testing estimation performance under conditions of ISS 
disturbance frequency uncertainty can be found in Figure 4-61. 
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Figure 4-61. Monte Carlo Results Testing ISS Disturbance Frequency 
Uncertainty (no spin up) 
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From Figure 4-58 it is evident that the Kalman filter performance is insensitive to rotor 
disturbance amplitude per the defined distributions. This result is important, because 
rodent motion and rodent mass modeling discrepancies are expected during normal 
operations due to the fact that there is no way to predetermine the rodent motion or to 
predict the rodent mass fluctuations over extended study periods. 

The Monte Carlo analysis result, considering a time-varying plant and disturbance inputs, 
can be found in Figure 4-59. 


% Amplitude Error (x . ) n ■ 0.861 I a m 0.007 % Amplitude Error (y _ ) |i" 0.216 / o ■ 0 

I Of Of I Of Of 




% Amplitude Error (♦ xrotof ) i*.- 0.996 / o - 0.165 % Amplitude Error ( 4 .^ 1 .463 / o - 0.505 




Figure 4-59. Monte Carlo Results Testing Rotor Disturbance Amplitude 
Uncertainty (with spin-up) 


Similar results are seen between the time-invariant and time-varying Monte Carlo analysis, 
supporting the conclusion that the Kalman Filter is effective for both cases. 
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a sinusoidal approximation of the effect of a pulse train through the system dynamics. 
Modeling the ISS disturbances in such a manner introduced another complexity (i.e. the 
modeling of a pulse train disturbance input) into the CR state estimation example. 
Observability for both the plant and filter models was evaluated. Since both models are 
time varying, observability needed to be checked for each variation of the model dynamics 
using both a time-invariant test as well as a time-varying test for observability. It was 
shown that the plant model is observable over the entire range of spin frequencies, and the 
filter model is observable over the entire range of the combinations of constant spin and 
ISS disturbance frequencies, except for low spin frequencies below 0.015 Hz. This is 
reasonable since low spin frequencies would create very small rotor disturbances, thus not 
affecting the system. Also, the discrete Kalman filter equations and algorithm are 
introduced along with a method for calculating initial Kalman filter parameters. 

In Chapter 4 the results from testing using the solution method proposed in Chapter 3 were 
presented. After determining performance measures, comparisons were made between 
different fidelities of disturbance modeling to show that it is necessary to model all of the 
disturbances which act on the plant in order to achieve good estimation. Also, comparisons 
were made between two filter models (ISS Model 1 and ISS Model 2), with different levels 
of ISS disturbance frequency modeling, show that expanding the disturbance model within 
the filter model will increase the range of disturbance frequencies for which the Kalman 
filter is effective. Finally Monte Carlo analysis shows that errors are sensitive to plant 
uncertainties as well as ISS disturbance frequency uncertainties. Therefore, care must bet 
taken to ensure that the filter model is as close to the plant model as possible. An on-line 
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system identification process may be necessary to detect and automatically adjust for any 
changes between the plant and the filter models [33][34][35]. Monte Carlo analysis also 
showed that it will be important to expand the filter model to cover the entire range of ISS 
disturbance frequencies. 

In conclusion, the use of disturbance modeling within the filter dynamics has proven to be 
useful in situations where the disturbance inputs into the plant are not available. Since 
many user defined parameters are not changed through the testing process, namely the 
process noise covariance, Q, and the initial error covariance, Po, the results provided are not 
the best possible. Therefore, for future work a method for updating Q and calculating Po in 
an optimal fashion should be investigated. 
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